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SUMHARY 


An  analysis  based  on  a  two-layer  interaction  model  is  presented  for 
the  prediction  of  transonic  shock-wave,  boundary- layer  interaction  with  em¬ 
phasis  on  turbulent  separated  flow.  In  this  analysis  finite  difference 
techniques  are  used  to  solve  the  viscous  layer  equations,  expressed  in  a 
defect  form,  and  the  stream  function-vorticity  representation  of  the 
inviscid  flow.  A  global  technique  is  used  to  iteratively  solve  these 
coupled  sets  of  equations.  Normal  pressure  gradients  and  imbedded  shock 
effects  are  included  in  this  analysis.  There  are  two  major  conclusions 
which  can  be  drawn  from  the  present  work.  First,  favorable  comparisons 
which  have  been  obtained  with  the  separated  data  of  Rooi  (M„  «  1 .4)  demon¬ 
strate  that  the  present  analysis  is  capable  of  accurately  resolving  many  of 
the  details  of  transonic  shock-wave,  boundary- layer  interaction.  Second, 
these  results  show  that  for  transonic  shock  induced  separation  the  effect 
of  displacement  thickness  interaction  dominates  over  that  produced  by 
imbedded  shock  effects  and  normal  pressure  gradients.  Calculations  made 
with  a  modified  algebraic  turbulence  model  demonstrate  that  for  separated 
cases  the  computed  results  are  more  sensitive  to  the  turbulence  model  than 
to  whether  or  not  normal  pressure  gradients  are  included.  Additional 
results  are  also  presented  for  transonic  shock-wave,  boundary- layer 
interaction  using  the  stream  function-vorticity  inviscid  technique  to 
implement  the  concepts  of  rotational  inviscid  flow  theory  for  the 
representation  of  transonic  interactions.  Finally,  a  new  three-layer 
interaction  model  is  presented  with  computed  results  shown  for  low  speed 
separated  flow  interactions. 
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INTRODUCTION 


The  prediction  of  the  flow  field  in  transonic  shock-wave,  boundary- 
layer  interaction  is  an  aerodynamic  problem  of  immense  practical  importance 
in  both  external  aerodynamics  and  turbomachinery  applications  since  this 
flow  phenomena  has  a  large  influence  on  the  overall  aerodynamic  forces  and 
often  determines  the  upper  limit  in  performance  of  each  aerodynamic  compo¬ 
nent,  Examples  of  such  components  include  airfoils  in  supercritical  flow- 
botb  external  and  internal  (turbomachinery  applications),  exhaust  nozzles, 
and  engine  inlets.  These  interactions  can  lead  to  premature  separation 
which  can  result  in  buffeting,  large  increases  in  drag,  and  reduced  engine 
performance.  In  addition,  it  is  difficult  to  predict  full  scale  shock 
interaction  phenomena  from  small  scale  wind  tunnel  data  because  of  the 
current  inability  to  match  experimentally  full  scale  Reynolds  number  flows. 
Hence,  there  is  a  strong  need  to  develop  analytical  and  computational  tech¬ 
niques  for  the  analysis  of  the  problem  of  transonic  shock  wave  boundary- 
layer  interaction. 

The  intersection  of  a  shock  wave  with  a  turbulent  boundary  layer  in 
transonic  flow  results  in  a  complex  flow  pattern.  Figure  1  shows  a 
schematic  diagram  of  the  flow  structure  which  occurs  for  both  attached  and 
separated  flow.  In  both  cases,  disturbances  created  by  the  impingement  of 
the  shock  wave  on  the  boundary  layer  propagate  upstream  in  the  subsonic 
portion  of  the  boundary  layer.  This  increase  in  thickness  creates  compres¬ 
sion  waves  which  intersect  with  the  shock  wave  and  cause  it  to  bend 
forward.  The  increase  in  the  pressure  at  the  wall  is  a  steep,  but  continu¬ 
ous  distribution  since  it  lies  below  the  sonic  line.  Away  from  the  wall 
discontinuities  may  occur  in  the  viscous  layer  if  the  shock  penetrates  this 
region.  If  the  upstream  Mach  number  is  not  too  high  the  flow  remains 
attached  with  the  flow  pattern  shown  on  the  left  side  of  Fig.  1.  However, 
an  increase  in  upstream  Mach  number  to  a  value  of  1 .3  or  greater  will 
result  in  flow  separation  with  a  more  complex  pattern  as  shown  on  the  right 
side  of  Fig.  1.  Generally  flow  separation  results  in  a  bifurcation  of  the 
shock  wave  into  a  front  and  rear  leg  near  the  wall  with  a  vortex  sheet 
emanating  from  the  triple  point.  A  separation  bubble  of  slow  recirculating 
flow  occurs  at  the  foot  of  the  shock  wave  in  the  lower  portion  of  the  vis¬ 
cous  layer.  In  both  the  attached  and  separated  flow  cases  the  upstream  in¬ 
fluence  of  the  shock  wave  is  typically  about  5  undisturbed  boundary- layer 
thicknesses;  the  downstream  extent  of  the  interaction  depends  on  the  shock 
strength.  The  short  streamwise  length  scale  of  the  interaction  region 
results  in  a  rapid  turning  of  the  flow  which  is  generally  thought  to  induce 
large  normal  pressure  gradients  in  the  viscous  layer.  The  overall  objec¬ 
tive  in  the  present  investigation  has  been  to  develop  s  reliable  and 
accurate  method  for  the  prediction  of  the  detailed  properties  of  the  flow 
which  occur  when  a  transonic  normal  shock  wave,  of  sufficient  strength  to 
cause  flow  separation,  impinges  on  a  turbulent  boundary  layer.  Only  the 
flow  in  the  immediate  vicinity  of  the  shock-wave,  boundary- layer 
interaction  has  been  considered  since  the  focus  of  the  present  work  was  to 
use  the  analysis  developed  herein  to  give  a  better  understanding  of  the 
relative  importance  of  imbedded  shock  effects  and  normal  pressure  gradients 
in  transonic  shock  induced  separation.  Before  the  details  of  this  analysis 
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are  given,  a  brief  discussion  is  presented  on  previous  approaches  which 
have  been  developed  for  this  problem. 

The  interaction  of  a  shock  wave  and  a  turbulent  boundary  layer  has 
been  the  subject  of  numerous  experimental  and  theoretical  investigations  in 
the  past  forty  years.  Only  a  brief  review  of  some  of  this  work  will  be 
given  here;  for  a  siore  detailed  account  of  past  work  the  reader  is  referred 
to  the  recent  review  articles  by  Adamson  and  Messiter  (Ref.  1),  Melnik 
(Ref.  2),  and  Lock  and  Firmin  (Ref.  3). 

Among  the  experimental  studies  is  the  work  of  Seddon  (Ref.  4),  Rooi 
(Refs.  5  and  6),  and  the  more  recent  work  of  Om  (Refs.  7  and  8)  in  which 
detailed  measurements  were  made  of  the  flow  properties  in  the  immediate 
vicinity  of  a  normal  shock  wave  impinging  on  a  flat  plate  turbulent  boun¬ 
dary  layer.  In  these  studies  the  shock  wave  was  of  sufficient  strength  to 
separate  the  turbulent  boundary  layer.  Seddon  and  Kooi's  experiments  were 
conducted  in  a  planar  configuration  whereas  Om's  measurements  were  made  in 
an  axisymmetric  tunnel  to  avoid  the  possible  three-dimensional  effects 
which  can  occur  in  planar  channel  experiments.  Detailed  data  were  obtained 
in  these  experimental  studies  for  the  pressure  distributions,  velocity  pro¬ 
files  and  associated  integral  properties,  and  the  structure  of  the  flow 
field.  In  all  of  these  flow  studies  a  lambda  shock  pattern  was  observed 
above  the  locally  separated  turbulent  boundary  layer.  An  important  issue 
in  the  turbulent  shock  boundary- layer  interaction  problem  is  whether  or  not 
the  shock  wave  penetrates  the  boundary  layer,  that  is,  do  the  compression 
waves  which  result  when  the  shock  strikes  the  rotational  flow  coalesce 
inside  or  outside  of  the  viscous  layer?  In  the  experimental  studies  of 
Rooi  on  transonic  shock  wave  boundary- layer  interaction,  it  was  observed 
that  at  a  free  stream  Mach  number  of  1.40  compression  wave  coalescence 
occurs  well  outside  the  boundary  layer  to  form  the  forward  leg  of  the 
oblique  shock.  At  a  higher  Mach  number  of  1.54  East  (Ref.  9)  observed  that 
the  "oblique  shock  was  formed  at  the  edge  or  probably  even  inside  the 
boundary  layer”.  Clearly  the  development  of  a  theoretical  method  which 
allows  for  imbedded  shock  effects  and  normal  pressure  gradients  would 
increase  the  current  ability  to  quantitatively  predict  the  influence  of 
Mach  and  Reynolds  number  on  transonic  interactions. 

Theoretical  approaches  to  the  problem  of  transonic  shock-wave, 
boundary- layer  interaction  range  from  numerical  solutions  of  the  Reynolds 
averaged  Navier-Stokes  equations  to  the  use  of  asymptotic  methods  which 
attempt  to  delineate  the  detailed  siathematical  structure  of  the  flow  field. 
Examples  of  the  former  approach  are  the  computations  made  by  Viegas  and 
Horstman  (Ref.  10)  and  Shea  (Ref.  11).  This  work  as  well  as  similar 
efforts  of  others  has  shown  that  this  approach  can  predict  the  general 
features  of  the  interaction  region  but  to  date  a  preferred  technique  has 
failed  to  emerge  from  these  numerical  studies.  Fart  of  the  difficulty  is 
that  with  a  full  numerical  approach  the  grid  must  be  carefully  chosen  in 
order  to  capture  the  flow  phenomena  which  occurs  at  different  length  scales 
in  the  interaction  problem.  It  is  well  recognised  that  if  the  Reynolds 
number  is  sufficiently  high,  as  it  generally  is  for  most  aerodynamic 
problems  of  interest,  the  flow  field  can  be  divided  into  viscous  and 
inviscid  regions  even  if  flow  separation  is  present.  Analytical  approaches 
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iriiich  split  the  flow  field  into  several  regions  will  clearly  provide  the 
best  overall  resolution  provided  that  the  proper  formulation  and  matching 
conditions  can  be  established  for  each  of  these  regions.  This  is  the  goal 
of  asymptotic  methods  which  are  described  next  for  the  application  to 
transonic  shock-wave,  boundary- layer  interaction. 


In  recent  years  large  Reynolds  number  asymptotic  expansions  have  been 
used  to  develop  local  solutions  for  weak  shock  waves  interacting  with  tur¬ 
bulent  boundary  layers.  Among  these  efforts  is  the  work  of  Melnik  and 
Grossman  (Ref.  12),  Adamson  and  Feo  (Ref.  13),  and  Messiter  (Ref.  14)  in 
which  they  have  shown  that,  in  the  absence  of  flow  separation,  the  interac¬ 
tion  of  a  normal  shock  wave  with  a  turbulent  boundary  layer  can  be 
represented  by  the  interaction  between  the  outer  part  of  the  turbulent 
boundary  layer,  treated  as  an  inviscid  rotational  flow,  and  the  incident 
shock  wave.  For  weak  shock  waves,  they  have  concluded  that  the  thin  vis¬ 
cous  sublayer  near  the  wall  plays  a  passive  role  in  the  overall  interaction 
process.  These  asymptotic  theories  have  evolved  from  the  pioneering  work 
of  Lighthill  (Ref.  15);  Melnik  (Ref.  2)  provides  an  excellent  discussion 
of  the  ensuing  work  between  Lighthill* s  paper  and  the  present-day  work  in 
asymptotic  analysis  of  transonic  interactions. 


Melnik  and  Grossman  (Ref.  12)  obtained  numerical  solutions  to  the 
inviscid  nonlinear  governing  equations  which  resulted  from  their  asymptotic 
analysis.  Their  solutions  demonstrate  that  the  mechanism  which  causes  the 
shock  to  bend  forward  and  blend  smoothly  into  the  sonic  line  is  contained 
entirely  in  the  interaction  between  the  inviscid  rotational  flow  and  exter¬ 
nally  imposed  shock  wave.  Later  Melnik  and  Grossman  (Ref.  16)  showed 
favorable  comparison  with  the  wall  pressure  determined  experimentally  by 
Gadd  (Ref.  17)  which  demonstrated  the  accuracy  of  their  asymptotic  method 
for  describing  the  flow  in  a  case  where  flow  separation  does  not  occur.  If 
flow  separation  does  occur,  which  is  the  subject  of  the  present  report, 
further  work  is  needed  in  the  development  of  asymptotic  methods  in  order  to 
provide  a  rational  theory  for  shock  induced  turbulent  separated  flow. 


The  approach  taken  in  the  present  investigation  is  an  outgrowth  of 
interacting  boundary  layer  theory  (IBLT)  which  is  now  briefly  reviewed.  In 
recent  years  significant  progress  has  been  made  in  the  development  of 
methods  for  the  iterative  solution  of  the  viscous  and  inviscid  equations 
through  displacement  thickness  coupling.  In  the  initial  development  period 
of  IBLT,  work  was  principally  focused  on  establishing  inverse  boundary 
layer  methods,  beginning  with  the  work  of  Catberall  and  Mangier  (Ref.  18) 
to  eliminate  the  classical  singularity  that  occurs  at  the  separation  point 
when  the  boundary- layer  equations  are  solved  subject  to  a  prescribed 
pressure  distribution.  Numerous  inverse  boundary- layer  methods,  discussed 
in  Appendix  A,  have  been  developed  since  this  earlier  work.  Following  this 
period,  work  was  directed  toward  the  development  of  procedures  which  permit 
inverse  boundary- layer  and  inviscid  solution  techniques  to  be  iteratively 
linked  together.  Two  methods  which  have  received  wide  use  are  the 
semi  -inverse  schemes  of  LeBalleur  (Ref.  19)  and  Carter  (Ref.  20).  These 
schemes  have  been  successfully  used  in  a  number  of  applications  of  IBLT  for 
the  prediction  of  laminar,  transitional,  and  turbulent  separated  flow. 


Current  efforts  in  the  development  of  IBLT  began  several  years  ago 
with  the  work  of  LeBalleur  (Ref.  21)  in  which  he  described  a  generalized 
viscous  layer  approach  where  the  governing  equations  are  expressed  as  the 
difference  (or  defect)  between  the  viscous  and  inviscid  flow  equations. 

This  treatment  permits  the  inclusion  of  normal  pressure  gradients  and 
embedded  shock  effects  in  the  solution  of  shock-wave,  boundary- layer 
interactions.  LeBalleur  (Ref.  21),  East  (Ref.  22),  and  Lock  and  Firmin 
(Ref.  3)  have  utilized  these  concepts  to  develop  generalized  integral 
boundary  layer  approaches  which  include  these  effects. 

The  present  work  is  based  on  a  viscous-inviscid  interaction  flow  model 
which  is  described  in  detail  in  the  next  section  of  this  report.  This 
model  principally  consists  of  a  differential  treatment  of  the  defect 
formulation  laid  out  by  LeBalleur  (Ref.  21)  utilizing  a  generalization  of 
the  inverse  boundary- layer  method  of  Carter  (Ref.  23)  coupled  to  the  stream 
function-vorticity  inviscid  method  of  Hafez  and  Lovell  (Ref.  24). 
Application  of  this  generalized  IBLT  formulation  to  several  problems  of 
transonic  shock-wave,  boundary- layer  interaction  are  presented  in  the 
section  on  Results  and  Discussion.  The  experimental  data  of  Kooi  (Ref.  5) 
is  used  to  evaluate  this  approach.  Numerous  conclusions  are  dra%m  from  the 
present  study;  the  principal  one  is  that  displacement  thickness  coupling  is 
the  dominant  effect  which  must  be  included  in  the  prediction  of  separated 
transonic  shock-wave,  boundary- layer  interaction. 


SYMBOLS 


Coefficient  of  specific  heat  at  constant  volume 
Viscous  perturbation  stream  function 
Velocity  ratio 
Total  enthalpy  ratio 

Displacement  thickness  divided  by  boundary  layer  thickness 
Total  enthalpy 

Approximate  viscous  layer  streamline  curvature 
Perturbation  mass  flow 
Mach  number 
Static  pressure 

Difference  between  pressure  in  real  viscous  flow  and  equivalent 
inviscid  flow 

Total  pressure 

Prandtl  number 

Reynolds  number 

Gas  constant 

Entropy 

Static  temperature 

Static  temperature  ratio 

Velocity  component  in  x-direction 

Velocity  component  in  y-direction 

Coordinate  measured  along  flat  plate 

Coordinate  measured  perpendicular  to  flat  plate 

Inviscid  velocity  gradient  parameter 

Ratio  of  coefficients  of  specific  heat 


Density 
Density  ratio 
Stream  function 

Inviscid  perturbation  stream  function 
Vorticity 

Relaxation  factor  in  global  iteration  scheme 

Molecular  viscosity  coefficient 

Eddy  viscosity  coefficient 

Boundary  layer  thickness 

Momentum  thickness 

Displacement  thickness 

Transformed  x-coordinate 

y'-Coordinate  divided  by  6* 

y-Coordinate  divided  by  6 

Free  stream  conditions 

Conditions  at  upstream  boundary 

Inviscid  quantity  at  reference  surface 

Inviscid  quantity  at  the  displacement  surface 

Inviscid  quantity  at  the  wall 

Quantity  in  equivalent  inviscid  flow 

Edge  of  boundary  layer 

Inviscid  prediction  in  global  iteration  technique 
Viscous  prediction  in  global  iteration  technique 
Incompressible 


VISCOUS-INVISCID  INTERACTION  FLOW  MODEL 


In  contrast  with  the  asymptotic  theory  developed  by  Melnik  and  Gross- 
man  (Ref.  12)  and  others  for  shock-wave,  turbulent  boundary- layer  inter¬ 
action  for  attached  transonic  flow,  a  rational  theory  has  yet  to  be 
developed  for  the  strong  interaction  case  in  which  the  shock  wave  is  of 
sufficient  strength  to  separate  the  boundary  layer.  In  the  absence  of  such 
a  formal  theory,  the  present  study  has  been  directed  toward  the  development 
of  two  viscous-inviscid  interaction  flow  models  for  strongly  interacting 
flows  with  separation.  These  models  are  considered  reasonable  since  their 
construction  is  a  natural  outgrowth  of  the  knowledge  which  does  exist  on 
the  nultiple  layer  structure  of  laminar  and  turbulent  interacting  flows. 

These  models  are  shown  schematically  in  Fig.  2.  In  both  models  the 
Reynolds  number  is  assumed  to  be  large  and  therefore  the  viscous  effects 
are  confined  to  a  relatively  thin  layer  near  the  surface.  The  main 
difference  in  these  two  models  is  the  treatment  of  the  viscous  layer.  In 
the  two-layer  model  the  boundary- layer  equations,  including  an  approximate 
treatment  of  normal  pressure  gradients,  are  solved  across  the  entire 
viscous  layer;  exterior  to  this  region  the  inviscid  equations  are  solved. 

In  this  model  the  inviscid  flow  may  be  rotational  due  to  an  upstream 
nonuniformity  or  the  vorticity  produced  by  a  curved  shock.  Matching 
conditions  are  imposed  which  result  in  a  smooth  merging  of  the  viscous  and 
inviscid  layers.  In  the  absence  of  normal  pressure  gradients  this  model 
reduces  to  that  of  classical  displacement  thickness  interaction.  The 
three  layer  model  consists  of  an  outer  inviscid  region,  which  is  coupled  to 
a  rotational  inviscid  layer  to  represent  the  outer  portion  of  the  upstream 
turbulent  profile,  and  a  viscous  sublayer  near  the  wall  to  accommodate  the 
imposition  of  the  no  slip  conditions.  In  this  model,  as  in  the  two-layer 
model,  rotational  effects  may  occur  in  the  outer  inviscid  region;  however 
in  this  model  there  is  an  additional  thin  vortical  layer  near  the  surface 
whose  thickness  is  of  the  order  of  the  upstream  boundary-layer  thickness. 

In  the  Melnik-Grossman  (Ref.  12)  asymptotic  analysis  for  attached  flows  the 
viscous  sublayer  plays  a  passive  role  in  the  overall  interaction  process 
and  thus  does  not  have  to  be  computed  unless  the  skin  friction  distribution 
is  desired.  For  separated  flows,  this  viscous  sublayer  plays  an  active 
role  and  therefore  its  effect  must  be  included  in  order  to  properly 
displace  the  outer  vortical  layer. 

The  bulk  of  this  report  is  directed  toward  the  development  and  verifi¬ 
cation,  through  experimental  comparisons,  of  the  two-layer  model.  The 
three-layer  model  is  presented  in  Appendix  A  along  with  computed  examples 
for  incompressible  separated  flow  which  demonstrate  this  method.  In  both 
of  these  models  the  inviscid  flow  has  to  be  capable  of  treating  rotational 
flow  due  to  an  upstream  nonuniformity,  curved  shock  wave,  or  the  vorticity 
contained  in  the  outer  portion  of  an  upstream  turbulent  boundary  layer.  In 
the  next  section  of  this  report  an  inviscid  analysis  based  on  the  stream 
function-vorticity  approach,  is  presented  for  the  computation  of  rotational 
transonic  flow.  Next  the  viscous  region  analysis,  used  in  the  two-layer 
model,  is  presented  followed  by  a  discussion  of  the  viscous-inviscid 


interaction  procedure  which  is  used  to  couple  the  inviscid  and  viscous 
flow  analyses. 


Inviscid  Analysis 

The  governing  equations,  boundary  conditions,  and  numerical  technique 
of  the  inviscid  analysis  used  in  the  present  study  on  transonic  shock  wave 
boundary  layer  interaction  are  described  here.  This  method  is  based  on  the 
compressible  stream  function  method  of  Hafez  and  Lovell  (Ref.  24)  which 
they  have  developed  for  the  analysis  of  transonic  inviscid  flow.  Prior  to 
their  work  difficulties  had  been  encountered  in  the  use  of  the  stream 
function  approach  in  the  transonic  speed  range  since  the  density  is  a 
double  valued  function  of  the  mass  flux  (there  is  a  subsonic  and  supersonic 
branch).  Hafez  and  Lovell  presented  a  method  which  overcomes  this 
difficulty  thereby  allowing  solutions  to  the  stream  function  equation  to  be 
obtained  for  transonic  flows.  This  approach  has  a  major  advantage  over 
that  based  on  the  potential  flow  approximation  since  it  is  equivalent  to  a 
solution  of  the  Euler  equations  provided  the  entropy  and  vorticity  produced 
by  a  curved  shock  or  an  upstream  nonuniformity  are  taken  into  account.  A 
second  advantage  of  the  stream  function  approach  is  that,  in  conservation 
form,  it  closely  resembles  the  form  of  the  full  potential  equation  for 
which  numerous  computational  techniques  have  been  developed  during  the  past 
ten  years. 

Governinx  Eauations  and  Boundarv  Conditions 


Since  the  present  work  is  focused  on  the  problem  of  a  two-dimensional 
normal  shock  wave  impinging  on  a  flat  plate  turbulent  boundary  layer,  the 
presentation  of  the  inviscid  equations  will  be  made  in  terms  of  the 
Cartesian  coordinates  x  and  y,  and  the  respective  velocities,  u  and  v. 
Introduction  of  a  stream  function  4'  defined  by 

dij/  dxl/ 

automatically  satisfies  the  continuity  equation 


dpu  dpv 

“7 —  +  "7 —  =  0 

dx  dy 

A  second-order  equation  is  determined  for  the  stream  function  by  the 
substitution  of  Eq.  (1)  into  the  definition  of  vorticity  u  given  by 

du  dv 
"  dy  ”  dx 


which  leads  to 
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I  d\py 


The  ioviscid  foraulation  is  completed  by  the  definition  of  the  total 
enthalpy,  assumed  constant  in  the  present  study. 


the  equation  of  state  written  in  the  form 


and  the  Crocco  relation 


P  ' 

Equation  (7)  shows  that  the  ratio  of  vorticity  to  pressure  is  constant 
along  a  streamline  since  s  >>  s(<|^)  only. 

The  definitions  given  in  Eq.  (1)  and  Eqs.  (3)'-(7)  constitute  seven 
equations  for  the  seven  unknowns  ip,  u,  s,  p,  p,  u,  and  v.  These  equations 
have  been  expressed  in  nondimen sional  form  in  which  the  x-  and  y-coordinqtes 
have  been  nondimensionalized  by  the  incoming  boundary- layer  thickness,  6^ 
shown  schematically  in  Fig.  3,  the  velocities  by  the  free  stream  value, 
u„,  the  density  by  Cs-i  and  the  pressure  by  PooUoo  where  an  asterisk  is  used 
here  to  denote  dimensional  quantities. 


The  boundary  conditions  imposed  for  the  solution  of  these  equations 
are  shown  schematically  in  Fig.  3  for  the  analysis  of  a  normal  shock  wave 
impinging  on  a  flat  plate  boundary  layer.  At  the  upstream  boundary 
all  dependent  variables  are  prescribed  from  the  incoming  free  stream 
flow  with  a  rotational  layer  near  the  wall  to  represent  the  outer  portion 
of  the  turbulent  boundary  layer.  At  the  outer  boundary  a  normal  shock  wave 
is  imposed  by  setting  the  values  of  p  and  u  from  the  Rankine-Hugoniot 
normal  shock  wave  conditions.  In  the  "free  air"  calculations,  with  the 
outer  boundary  placed  sufficiently  far  from  the  wall  so  that  the  flow  is 
undisturbed  at  y  *  y  ,  the  values  of  p,  u,  and  4'  are  held  constant 
along  this  boundary .°‘^{n  calculations  which  were  made  to  compare 
the  present  method  with  experimental  data,  it  was  found  necessary  to 
impose  conditions  along  this  outer  boundary  to  represent  the  blockage 
effect  which  occurs  when  a  shock-wave, boundary -layer  interaction  occurs  in 


a  channel.  In  this  case  the  outer  boundary  is  not  a  streamline  but  rather 
the  condition 


(8) 

is  imposed  where  pu  is  prescribed  from  the  experimental  measurements. 

At  the  wall  the  presence  of  the  viscous  flow  is  represented  in  the 
inviscid  analysis  as  a  Dirichlet  condition  on  the  stream  function,  <l^(x,0)  > 
il<v(x).  The  relationship  between  '{<v(x)  and  the  viscous  solution  will  be 
discussed  in  detail  in  the  section  on  the  Viscous-Inviscid  Interaction 
Procedure.  At  the  do%mstream  or  outflow  boundary  is  set  to  zero; 
calculations  were  also  made  with  i|(  *  0  at  this  boundary  with  only  small 
changes  in  the  results. 

Numerical  Technique 

In  this  section  the  numerical  technique  used  to  solve  the  inviscid 
formulation  is  presented.  This  technique  closely  follows  that  presented 
previously  by  Hafez  and  Lovell  (Ref.  24)  for  the  computation  of  transonic 
flow  with  the  stream  function-vorticity  approach. 

Coordinate  stretchings,  X  -  g(x)  and  Y  >  f(y),  have  been  introduced  to 
place  fine  grids  near  the  shock  wave  and  the  wall  to  resolve  the  large 
gradients  which  occur  in  these  regions.  The  transformed  stream  function- 
vorticity  equation  is 


^  ax  (  n  ax  /  ^  dY  \p  av/  ■ 
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(9) 


where  a  prime  is  used  to  denote  differentiation.  Equation  (9)  is  represen¬ 
ted  with  central  difference  approximations  which  results  in 
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for  the  computational  molecule  shown  in  Fig.  4.  With  o  and  u  assumed 
known,  Eq.  (10)  is  solved  with  the  zebroid  scheme  (Ref.  25)  which  utilizes 
line  relaxation  along  alternating  black  and  white  horizontal  lines.  A 
finite  difference  approximation  to  a  term  has  been  explicitly  added  to 
this  equation  to  guarantee  convergence  in  the  supersonic  zone  (Ref.  26). 


After  the  stream  function  is  computed  the  u-component  of  velocity  is 
deduced  by  the  integration  of  Eq.  (3)  from  the  outer  boundary  Y  >  1  where  u 
is  prescribed.  This  results  in 

Y 

u(X,Y)  =  u(X,l)-  -Jy  /  (-cu  +  g'  dY  (11) 

^  Yal 


Vitb  b)  and  v  aasumed  known,  Eq.  (11)  ia  integrated  with  the  trapezoidal 
rule  to  yield  u  at  the  segment  centers  denoted  by  an  x  in  Fig.  4.  The  v- 
component  of  velocity  is  then  calculated  at  the  cell  centers,  denoted  by 
a  ■»’  in  Fig.  4,  with  the  finite  difference  representation  of 


g'ta»///dx) 

V  =  -  -  u 

f  (d\l//dr) 


(12) 


In  Eq.  (12)  the  value  of  u  at  the  cell  center  is  obtained  by  averaging  the 
segment  center  values. 


The  density  is  then  updated  at  the  cell  centers  using  a  combined  form 
of  the  total  enthalpy  and  state  equations  which  gives 


I 
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(13) 


The  values  at  the  segment  centers  are  computed  by  averaging  the  cell  center 
values.  It  is  well  known  that  in  order  to  maintain  numerical  stability  in 
transonic  flow  computations  it  is  necessary  to  add  an  artificial  viscosity 
term  in  the  supersonic  region.  In  the  present  calculations  the  retarded 
density  method  of  Hafez,  et  al.  (Ref.  26)  has  been  used  to  add  artificial 
viscosity  using  a  modified  density 


AV 

with 

/I  =  max  [o,c  ~  ^  ] 

where  c  is  an  assigned  constant  and  M  is  the  local  Mach  number.  The 
term  is  computed  with  an  upwind  finite  difference  approximation.  The  use 
of  artificial  viscosity  in  the  inviscid  analysis  causes  the  shock  wave  to 
be  spread  over  several  grid  points;  calculations  will  be  presented  in  the 
Results  and  Discussion  section  which  show  that  this  shock  smearing  is 
significantly  smaller  than  that  which  results  due  to  the  real  viscous 
effects. 


(14) 


(15) 
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The  step*  in  this  numerical  procedure  are  completed  by  using  the 
Crocco  relation,  Eq.  (7),  to  update  the  vorticity  after  every  10  sweeps  of 
the  inviscid  flow  equation.  In  the  present  calculations  the  entropy  rise 
at  the  shock  has  been  neglected;  hence  the  entropy  and  vorticity 
distribution  is  due  solely  to  that  specified  at  the  upstream  boundary  to 
represent  the  outer  portion  of  the  incoming  turbulent  boundary  layer. 

Since  this  distribution  is  known  a  priori  an  analytical  curve  fit  is  used 
to  represent  s  ■  *(<!').  In  the  numerical  calculations,  with  the  latest 
values  of  ^  and  P  assumed  known  at  each  grid  point,  s  and  ds/di|»  are 
computed  analytically  from  this  curve  fit  and  the  value  of  u  is  updated  at 
each  point  by  the  substitution  of  Eq.  (6)  for  the  pressure  into  the  Crocco 
relation,  Eq.  (7). 

Nelgect  of  the  entropy  rise  at  the  shock  wave  does  not  result  in  a 
serious  error  in  the  present  calculation  since  the  Rankine-Hugoniot 
relations  have  been  used  to  specify  the  correct  shock  jump  conditions  for  p 
and  u  at  the  outer  computational  boundary.  Free  air  calculations  at  Mo,  ■ 
1.4  have  shown  that  the  value  of  the  pressure  behind  the  normal  shock  wave 
is  overpredicted  by  only  about  5  percent.  However,  as  will  be  discussed  in 
greater  detail  in  the  Results  and  Discussion  section,  this  approximation 
does  prevent  the  formation  of  a  lambda  shock  in  the  calculation  which 
generally  occurs  when  the  shock  wave  is  of  sufficient  strength  to  separate 
the  boundary  layer. 


Viscous  Analysis 

In  this  section  the  formulation  and  solution  technique  is  given  for 
the  viscous  region,  which  is  contained  in  the  two-layer  viscous-inviscid 
interaction  flow  siodel  discussed  previously.  This  analysis  is  an  extension 
of  the  first-order  inverse  boundary- layer  analysis  (Ref.  23)  to  now  include 
normal  pressure  gradients  and  asymptotic  merging  with  the  outer  inviscid 
flow  at  the  viscous  layer  edge. 

Covernine  Equations  and  Boundary  Conditions 

The  governing  equations  for  the  viscous  flow  analysis  are  written  as 
follows  in  terms  of  Cartesian  coordinates  and  the  corresponding  velocity 
components,  u  and  v,  as  shown  in  Fig.  5: 


(18) 


dp  =  -  2 
T"  'f/’U 
dy 
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Equations  (16)(  (17)  and  (19)  are  the  usual  continuity  equations  and  boun¬ 
dary-layer  foms  of  the  x-nomentum  and  total  enthalpy  equations,  respec¬ 
tively*  Equation  (18)  is  an  approximate  form  of  the  normal  momentum 
equation,  suggested  by  LeBalleur  (Ref.  21)  and  Lock  and  Firmin  (Ref.  3),  in 
which  it  has  been  assumed  that  the  viscot^s  terms  are_negligible  and  the 
convection  terms  can  be  approximated  by  <Pu^  where  <(x)  is  an  average 
viscous  layer  curvature. 

Additional  equations  are  needed  to  complete  the  specification  of  the 
viscous  layer  formulation.  The  molecular  viscosity  coefficient  is  given  by 
the  Sutherland  relation  which  for  air  is  written  in  nondimensional  form  as 


t98.6“R 
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where  T,  and  denote  the  dimensional  (expressed  in  degrees  Rankine)  and 
nondimensional  reference  temperatures  of  the  free  stream,  respectively. 
The  equation  of  state  is  written  as 

p  -  (21) 


where  the  temperature  is  deduced  from  the  total  enthalpy  definition 


H  =  T  + 


Equations  (16)-(22)  constitute  a  set  of  seven  equations  in  the  seven 
unknowns,  u,  v,  p,  P,  H  and  T,  provided  that  the  average  curvature  is 
known  and  that  the  Reynolds  stresses  are  related  to  the  mean  flow  through  a 
turbulence  model.  In  this  investigation  the  Reynolds  stresses  have  been 
related  to  the  mean  flow  by  an  eddy  viscosity  representation  given  by 
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(24) 
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In  the  present  formulation  a  constant  turbulent  Prandtl  number  Prt is 
assumed,  which  is  defined  in  the  same  manner  as  Cebeci  and  Smith  (Ref.  27) 


thereby  giving 
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The  turbulence  models  used  in  the  present  calculations  will  be  discussed 
later. 

The  variables  in  Eqs.  (16)-(26)  are  expressed  in  nondimensional  form 
with  lengths  scaled  by  the  upstream  boundary- layer  thickness  6^  as  shown  in 
Fig.  5,  velocities  by  the  free  stream  value,  u^,  pressure  by  twice  the  free 
stream  dynamic  pressure,  o^u*^,  temperature  by  u„  /cp,  total  enthalpy  by 
u*^,  density  by  o*,  and  viscosity  coefficients  (u  and  e:)  by  u*.  An 
asterisk  is  used  here  and  in  Fig.  5  to  denote  dimensional  quantities.  The 
Reynolds  number  which  appears  in  Eqs.  (17)  and  (19)  is 


Re« 


(27) 


The  equations  presented  above  describe  the  flow  in  the  viscous  layer 
adjacent  to  the  wall,  y  >  0.  As  stated  previously,  it  is  assumed  that  the 
Reynolds  number  is  sufficiently  high  so  that  the  overall  flow  field  can  be 
divided  into  a  viscous  layer  near  the  surface  and  an  inviscid  region 
outside  of  the  viscous  layer.  Hence,  at  a  sufficiently  large  distance  from 
the  surface,  the  viscous  solution  is  required  to  merge  smoothly  into  that 
for  the  inviscid  flow.  Therefore,  it  is  convenient  to  write  the  viscous 
layer  equations  as  the  difference  between  the  real  viscous  flow  (RVF)  and 
the  equivalent  inviscid  flow  (EIF)  where  the  latter  results  when  the 
inviscid  flow  is  solved  in  the  interior  of  the  viscous  region.  LeBalleur 
(Ref.  21)  pioneered  this  idea  of  writing  the  viscous  equations  in  this 
defect  form;  Lock  and  Firmin  (Ref.  3)  have  further  explored  this  idea  and 
introduced  the  concepts  of  RVF  and  EIF  mentioned  above. 


The  defect  form  of  the  viscous- layer  equations  used  in  the  present 
study  are  written  as  follows  in  which  the  i-subscript  is  used  to  denote  the 
variables  in  the  EIF: 
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It  has  been  found  convenient  in  the  viscous  layer  analysis  to  replace  the 
v-component  of  velocity  with  the  stream  function  using  the  relation 


=  -pv 


with  the  continuity  equation,  Eq.  (16),  now  replaced  by  Eq.  (28).  The 
difference  between  the  viscous  and  inviscid  pressures  is  denoted  by  0  ■ 
p-p^  in  Eqs.  (29)  and  (30).  It  is  observed  that  the  EIF  in  Eqs.  (28)-(31) 
has  been  expressed  in  terms  of  the  Euler  equations  instead  of  the  conser¬ 
vative  form  of  the  stream  function-vorticity  equation  discussed  in  the 
section  of  the  inviscid  flow  analysis.  This  step  is  taken,  including  the 
use  of  the  approximate  form  of  the  inviscid  normal  momentum  equation,  to 
insure  that  as  the  edge  of  the  viscous  layer  is  approached  the  imposition 
of  the  boundary  conditions 

u-ui  H-Hi  p-0  (33) 

results  in  a  smooth  merging  of  the  viscous  and  inviscid  solutions  as  shown 
in  Fig.  6  for  the  u-component  of  velocity. 


At  tbe  surface  y  ~  0  the  usual  viscous  boundary  conditions 


u='j'  =  0  H=H^(x)  or  -r—  -  QwM  (34) 

dy 

are  imposed,  where  denotes  the  nondimen sional  form  of  the  specified 

surface  heat  transfer. 

Transformed  Viscous  Formulation 

Based  on  previous  work  in  the  analysis  is  separated  boundary- layer 
flows,  the  finite  difference  solution  of  the  governing  equations  is 
facilitated  by  the  use  of  a  modified  form  of  the  Levy-Lees  transformation 
which  is  given  by: 

X 

^  ■/  ^ir  “ir «  =  JL  h  (  x )  = -r^h(x)  (3b) 

0  S 

In  Eq.  (35)  the  ir-subscripted  quantities  denote  the  inviscid  solution  at  a 
reference  surface  such  as  the  wall  or  the  displacement  body.  These  quanti¬ 
ties  depend  only  on  x;  further  discussion  on  their  treatment  in  the  present 
study  will  be  presented  in  the  section  on  the  Equivalent  Inviscid  Flow. 
Equation  (35)  shows  that  the  normal  coordinate  has  been  scaled  by  the  dis¬ 
placement  thickness.  Numerical  difficulties  which  were  encountered  in 
using  this  transformation  in  the  present  computations  of  transonic  shock 
induced  separated  flow  were  overcome  by  rescaling  the  normal  coordinate  by 
the  boundary- layer  thickness,  6  and  therefore  in  Eq.  (35),  h(x)  *  6  /i . 
Further  discussion  on  this  point  will  also  be  presented  in  the  section  on 
the  Equivalent  Inviscid  Flow. 

In  the  present  analysis  the  inviscid  flow  is  permitted  to  vary  over 
the  viscous  layer  thickness  due  to  the  strong  viscous-inviscid  interaction 
which  occurs  when  a  transonic  shock  wave  impinges  on  a  turbulent  boundary 
layer.  It  is  therefore  necessary  to  account  for  this  inviscid  variation  by 
generalizing  tbe  displacement  definition  to 

f  p.u.dy  =  Jj  pudy  (36) 

8'**^  o 

which  states  that  the  mass  flow  in  the  equivalent  inviscid  flow  is  equal  to 
that  in  the  real  viscous  flow.  In  the  present  investigation  it  has  been 
found  convenient  to  use  an  alternative  displacement  thickness  definition 
given  by 


- 


OO 

J  {pu^  -  pu)dy 


(37) 


Lock  (Ref.  28)  bee  explored  the  difference  in  these  two  definitions  of  6* 
and  has  concluded  that  they  differ  less  from  each  other  than  they  do  from 


the  classical  definition  obtained  by  setting  p.u. 
(37)  which  leads  to  ^  ^ 


in  Eqs. 


(36)  and 


(38) 


where  PgU^  is  constant  across  the  boundary  layer. 

In  the  viscous  layer  it  is  convenient  to  use  a  perturbation  stream 
function  to  represent  the  difference  between  the  viscous  and  inviscid 
stream  functions  written  as 


As  y  f  1  so  that  the  viscous  stream  function  merges  smoothly  with 
that  of  the  inviscid  solution.  At  the  surface,  y  -  0,  the  usual  boundary 
condition  *  0  is  imposed  by  setting  f  ~  0  since  at  the  wall  the  inviscid 
stream  function  is  set  to 


(x,o)=-m 


where  m  is  the  perturbation  mass  flow  given  by 


m:  p.  u-  S* 

ir 


(40) 


(41) 


The  boundary  condition  given  in  Eq.  (40),  which  is  the  representation  of 
the  VISCOUS  effects  in  the  inviscid  flow,  will  be  derived  later  in  the 
section  on  the  Viscous-Inviscid  Interaction  Procedure.  In  first  order 
theory,  where  the  inviscid  flow  velocity  is  invariant  on  the  scale  of  the 
boundary  layer,  the  inviscid  stream  function  varies  linearly  in  the  viscous 
region.  Hence  a  second  perturbation  stream  function,  is  introduced 
which  denotes  the  deviation  of  the  inviscid  stream  function  from  its  first 
order  linear  variation.  This  function  is  given  by 


-  m  (  ^-i) 


(42) 


which  when  combined  with  Eq.  (39)  yields  the  following  form  for  the  viscous 
stream  function 

»^rm(  TJ-f')  +  Ij, 

A  (43) 

Figure  7  shows  the  two  perturbation  stream  functions,  and  f,  and  the 
merging  of  the  inviscid  and  viscous  stream  functions  at  the  edge  of  the 
viscous  layer. 
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The  normalized  state  equation  becomes 


and  the  temperature  is  deduced  from  the  total  enthalpy  definition,  Eq. 
(22),  which  is  now  written  as 


The  v'component  of  velocity  is  deduced  from  the  transformation  of  Eq.  (32) 
given  by 


Ujj.  p  I- L  \h  j  i  h'h  S*/^  ^  (51) 


The  Mach  number  is  given  by 


(y-l)Ti, 


The  boundary  conditions  which  are  imposed  on  the  viscous  formulation  are 
given  by 

dg 

■7=0  F  =  f=0.  g  =  or  —  =  9^(0  (53^ 

OTj 


•7—00 


F  -F. 


g-Qi 


p-O 


where  g^(C)  is  a  prescribed  wall  heat  flux  distribution.  An  important 
feature  of  this  formulation  is  that  as  the  boundary- layer  edge  is 
approached  and  the  viscous  shear  and  heat  conduction  terms  vanish,  the 
viscous  flow  solution  is  required  to  merge  asymptotically  with  the  inviscid 
solution.  He^ce  the  boundary  conditions  given  in  Eq.  (54)  can  be  imposed 
at  any  value  n  >  1.  Numerical  tests  were  conducted  on  a  model  problem  and 
reported  earlier  in  Ref.  29  to  verify  this  asymptotic  behavior. 
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Equations  (45}-(30)  along  with  the  Sutherland  molecular  viscosity  law 
constitute  seven  equations  in  the  seven  unknowns  F,  £,  p,  g,  p,  T,  and  p, 
with  the  inviscid  solution  for  the  equivalent  inviscid  flow  assumed  known. 
The  eddy  viscosity  coefficient  c  is  related  to  these  flow  variables  through 
the  use  of  the  Cebeci-Smith  (Ref.  27)  algebraic  turbuence  model.  It  was 
found,  as  will  be  presented  in  the  section  on  the  RESULTS  and  DISCUSSION 
that  a  modification  to  the  outer  length  scale  in  this  model  improved  the 
agreement  of  the  predicted  results  with  transonic  flow  measurements. 

The  imposition  of  an  outer  boundary  condition  on  the  stream  function 
permits  these  equations  to  be  solved  in  the  inverse  form,  that  is,  the 
unknown  pressure  gradient  parameter  S  *  (1/u.  )/(du.  /d^)  is  deduced  along 
with  the  viscous  solution  for  a  prescribed  diitribuHon  of  the  perturbation 
mass  flow,  m.  This  approach  allows  the  present  generalized  viscous 
formulation  to  be  incorporated  into  the  interaction  iteration  procedure 
which  has  been  successfully  used  in  numerous  studies  (Ref.  30)  on 
first-order  viscous-inviscid  interaction.  It  is  observed  that  the  present 
analysis  is  easily  reverted  back  to  a  first-order  formulation  in  which  the 
pressure  is  constant  across  the  boundary  layer  by  setting  F.  =  p.  =  T.  = 

1,  *  0,  and  ^  »  0.  In  this  case  the  equivalent  invisid  ftow  is  ^ 

invariant  on  the  scale  of  the  boundary  layer  and  hence  u.~u,p.»p,p. 
■  Pg  where  the  edge  quantities,  denoted  by  the  e-subscript,  are  constant  a^ 
each  x-station  along  the  surface. 

In  the  present  study  several  approximations  have  been  made  to  simplify 
Eqs.  (45)-(50).  As  a  first  step  in  the  inclusion  of  normal  pressure 
gradients  in  an  interacting  boundary- layer  analysis  the  viscous  pressure 
has  been  set  equal  to  that  of  the  inviscid  flow  thereby  resulting  in  ^  =  0 
with  Eq.  (47)  eliminated  from  the  formulation.  This  step  is  thought  to  be 
a  reasonable  approximation  since  the  variation  of  the  pressure  across  the 
viscous  layer,  if  it  occurs,  will  exist  primarily  in  the  outer  region  as 
sketched  in  Fig.  8.  In  this  region  normal  pressure  gradients  are  induced 
in  the  inviscid  flow  due  to  the  rapid  turning  which  principally  occurs  in 
the  present  problem  at  the  foot  of  the  shock  wave.  Hence  we  allow  this 
inviscid  variation  of  pressure  to  be  imposed  on  the  viscous  solution  by 
setting  p(x,y)  -  p.(x,y)  which,  as  pointed  out  by  LeBalleur  (Ref.  21),  is 
less  restrictive  t^an  the  first  boundary- layer  approximation,  p(x)  » 
p.(x,0).  The  implementation  of  the  present  treatment  of  the  normal 
pressure  gradient  will  be  discussed  in  detail  in  the  next  section  on  the 
Equivalent  Inviscid  Flow.  In  the  future,  this  generalized  or  higher-order 
viscous-inviscid  interaction  flow  model  should  be  further  developed  by 
solving  Eq.  (47)  simultaneously  with  the  other  viscous  layer  equations. 

This  step  should  be  straightforward  with  the  present  inverse  solution 
procedure  since  the  mean  streamline  curvature,  <,  if  it  is  set  equal  to 
that  of  the  displacement  surface  as  recommended  by  LeBalleur  (Ref.  21)  and 
Lock  and  Firmin  (Ref.  3),  will  be  known  on  each  streamwise  calculation  of 
these  equations. 

In  addition  to  this  approximate  treatment  of  the  normal  pressure 
gradient,  the  present  work  has  been  focused  on  adiabatic  flows  with  the 
Prandtl  number  Pr  set  equal  to  unity..  This  approximation  allows  us  to  set 
g  ■  g.  >1  and  therefore  eliminate  the  total  enthalpy  equation,  Eq.  (48) 
from  the  present  analysis.  This  approximation  should  have  a  minimum  effect 


on  the  computed  results  since  in  adiabatic  flows  the  total  enthalpy  ratio 
g  only  varies  by  about  3  percent  across  the  viscous  layer.  The  inclusion 
of  the  total  enthalpy  form  of  the  energy  equation  to  study  the  effect  of 
heat  transfer  on  flows  with  strong  viscous-inviscid  interaction  is  straight 
forward. 


Equivalent  Inviscid  Flow 


The  inviscid  flow  solution  in  the  interior  of  the  viscous  layer,  which 
is  referred  to  as  the  equivalent  inviscid  flow  (EIF),  is  required  for  two 
purposes  in  the  present  viscous  layer  formulation.  First  it  is  needed 
since  the  viscous  equations  have  been  written  as  the  difference  (defect 
form)  between  the  EIF  and  the  real  viscous  flow  (RVF),  and  second,  in  the 
present  calculations  the  viscous  streamwise  pressure  gradient  has  been 
approximated  as  the  inviscid  streamwise  pressure  gradient  expressed  in 
generalized  Levy-Lees  variables  as 


dp  ,  . 

67  ‘ 


dx 


(x.y)  =- 


ir 
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r  _  dF, 


mh 


d  r  ,  ^ 


^  d-rj 


(55) 


Eq.  (55)  has  been  used  to  approximate  the  y-variation  of  the  x-pressure 
gradient  only  for  y  >  or  n  >  1  •  Below  y  *  6*  this  equation  has  been 
replaced  with  the  assumption  that  the  pressure  gradient  is  now  independent 
of  y  and  equal  to  that  of  the  EIF  at  the  displacement  thickness.  Hence  for 

y  <  6^' 


dp  i  . 
-7—  (x)  = 

OX 


dp. 


dx 


(x.8*) 


(56) 


which  is  shown  schematically  in  Fig.  8.  Equation  (56)  should  be  a 
reasonable  approximation  for  y  <  6^'  since  it  is  known  that  the  normal 
pressure  gradient  in  the  RVF  is  zero  at  the  wall  as  shown  in  Fig.  8. 

In  order  to  accommodate  this  treatment  of  the  pressure  gradient,  the 
inviscid  reference  surface,  about  which  the  viscous  formulation  is 
normalized,  has  been  chosen  as  the  displacement  surface  at  n  ~  1 .  Hence 
the  ir-subscripted  quantities  have  been  taken  as  the  inviscid  solution  at  y 
>  6*.  In  the  global  viscous-inviscid  interaction  procedure,  just  before 
the  viscous  solution  is  updated,  the  required  input  from  the  EIF  is 
computed.  Using  the  current  value  of  6*  linear  interpolation  in  the  y- 
direction  is  performed  at  each  x-station  |o  determine  the  inviscid  u-compo- 
nent  of  velocity  and  the  density  at  y  >  6  .  Denoting  these  values  as  u^^  ~ 
u,*  and  p.  ■  p  *,  respectively,  we  now  compute  the  ^-distribution  of 


using  the  trapezoid  rule  of  integration.  Belov  the  displacement  thickness 
F.  *  c.  “  1  and  thus  va  ■  0  in  this  region  so  that  the  y-dependence  of  the 
sireamvise  pressure  gradient  is  eliminated.  Thus  in  this  region  Eqs.  (35) 
and  (561  yield 


dx 


Ms*  ^S* 


(59) 


which  is  the  same  as  that  which  would  be  imposed  across  the  entire  viscous 
layer  in  a  tirst-order  analysis. 


It  was  mentioned  in  the  previous  section  that  it  was  found  necessary 
to  rescale  the  y-coordinate  with  the  local  boundary- layer  thickness  instead 
of  the  displacement  thickness  which  was  successfully  used  in  previous 
first-order  analyses  of  viscous-inviscid  interaction  (Ref.  30).  This  modi¬ 
fication  was  needed  since  some  caution  was  required  as  to  the  lateral 
extent  to  which  the  viscous  layer  equations  could  be  solved.  This  results 
because  the  viscous  layer  equations  have  been  written  in  nonconservative 
form;  numerical  difficulties  were  encountered  when  the  viscous  solution  was 
extended  too  far  outside  of  the  viscous  layer  edge  due  to  the  nearly  dis¬ 
continuous  behavior  in  the  flow  field  variables  in  that  region.  Rescaling 
by  ^  eliminated  this  problem  and  no  numerical  difficulties  were  encountered 
with  the  outer  boundary  in  the  viscous  solution  located  at  y/6  *  n  ^1.4. 


The  viscous  layer  equations  with  the  imposed  EIF  are  solved  with  an 
implicit  finite-difference  which  is  first-order  accurate  in  the 
C-coordinate  and  second-order  accurate  in  the  ^coordinate.  The  finite 
difference  representation  of  the  EIF  terms  is  treated  identically  to  the 
corresponding  viscous  flow  convection  terms  thereby  insuring  that  these 
terms  are  in  balance  as  the  edge  of  the  viscous  layer  is  approached.  Other 
than  the  inclusion  of  the  EIF  terms  the  numerical  treatment  is  identical 
with  that  used  in  the  firrt-order  viscous  analysis  which  is  presented  in 
Ref.  23. 


Viscous-inviscid  Interaction  Procedure 

The  viscous  induced  displacement  thickness  is  represented  in  the 
inviscid  formulation  as  a  Dirichlet  condition  for  the  stream  function  as 


I  M  '-V  'LV  "-1  ' 


shown  in  Fig.  3.  This  boundary  condition  is  discussed  in  this  section 
along  with  the  global  iteration  procedure  used  to  numerically  link  the 
inviscid  stream  function  analysis  to  the  generalized  viscous  formulation. 

Inviscid  Surface  Boundary  Condition 

The  defect  form  of  the  continuity  equation  is  written  as 

aiftvffv)  (JO, 

dx  dy 

Integration  of  Eq.  (60)  from  y  *  0  to  y  =  «=  with  the  imposition  of  the 
boundary  conditions 


y=0  ;>v=0  (61) 


1/>U  — yO.Uj  (62) 

yields 

^  oc 

(x,0)  =  —  ^(/.,u,-^u)dy  ,j3) 

Introduction  of  the  stream  function  from  Eq.  (1)  where  the  i-subscript  is 
used  here  to  denote  the  inviscid  value >  leads  to 


OW.  fj  y 

This  equation  is  integrated  from  x  ~  to  x  to  give 

»//.  (x.O)  -fj/.  (Xq.O)  =  -m  (x)  +  m(xp) 


(64) 


(65) 


where  m  is  the  perturbation  mass  flow  previously  defined  as 

X 

m(x)  =  J(/c>jUj-/Du)dy  (66) 

0 

It  is  convenient  to  choose  the  constant  of  integration  as  (x  ,o)  =  -m(x^) 
which  leads  to  the  final  result 
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as  the  inviscid  surface  boundary  condition. 


>bal  Iteration  Technique 


The  iteration  procedure  shown  in  Fig.  9  is  used  to  couple  the  tran¬ 
sonic  stream  function  analysis  to  the  viscous  formulation  described  in  this 
report.  This  iteration  technique  is  a  generalization  of  a  procedure  (Ref. 
20)  which  has  been  successfully  used  in  numerous  computations  of  first- 
order  viscous-inviscid  interaction.  Since  the  viscous  equations  are  now 
written  in  a  defect  form  the  viscous  and  inviscid  equations  are  much  more 
strongly  coupled  than  in  a  first-order  computation  where  the  pressure  is 
assumed  constant  across  the  boundary  layer.  Nonetheless  the  simple  fixed 
point  iteration  formula  shown  in  Fig.  9  has  been  successfully  used  in  the 
present  study  to  obtain  converged  interaction  solutions.  In  Fig.  9,  the 
update  formula  is  shown  with  the  relaxation  factor,  u,  set  to  unity;  in 
general  it  has  been  found  necessary  in  the  present  calculations  to  use 
underrelaxation  and  thus 


m''*’'  =  m' 


In  Eq.  (68)  u^«  and  u^.>  denote  the  viscous  and  inviscid  predictions, 
respectively o1^  the  inviscid  velocity  along  the  displacement  surface.  In 
the  computed  results  to  be  presented  in  the  next  section,  the  quantity, 

Ujfc^  has  been  approximated  as  the  inviscid  velocity  at  the  wall,  uy^.  This 
Approximation  was  made  to  eliminate  minor  oscillations  which  resulted  when 
u^jir  was  used  in  the  update  expression.  These  oscillations  are  probably 
duetto  the  interpolation  which  is  required  in  the  inviscid  calculation  to 
find  u^^  since  the  grid  is  aligned  to  the  wall  and  not  to  the  displacement 
surface.^  This  approximation  has  a  relatively  small  effect  on  the  results 
since  the  inviscid  flow  changes  only  slightly  between  the  wall  and  the 
displacement  thickness. 


RESULTS  AND  DISCUSSION 


la  this  section  results  are  presented  using  the  previously  discussed 
two-layer  model  (IBLT)  for  both  attached  and  separated  transonic  normal 
shock-wave,  boundary- layer  interaction  including  the  effects  of  normal 
pressure  gradients.  In  addition  calculations  are  presented  for  these  same 
cases  utilizing  an  inviscid  rotational  flow  model  to  represent  the  outer 
portion  of  a  turbulent  boundary  layer  interacting  with  an  incident  normal 
shock  wave.  Comparisons  made  between  these  two  approaches  delineate  the 
effect  of  ignoring  the  viscous  sublayer  in  the  rotational  inviscid  flow 
analysis.  Assessment  is  made  of  these  results  through  comparisons  with 
previous  analyses  and  experimental  data.  A  brief  discussion  is  presented 
of  an  improved  turbulence  model  which  was  applied  during  this  investigation 
and  has  been  found  to  give  improved  agreement  of  theory  and  experiment  for 
some  quantities  for  transonic  shock  induced  separated  flow.  Finally, 
calculations  are  shown  from  the  present  IBLT  which  assess  the  relative 
importance  of  boundary- layer  displacement  effects  and  normal  pressure 
gradients  for  transonic  shock  induced  separated  flow. 


Attached  Flow  Calculations 

A  number  of  investigators  (Ref.  2)  have  demonstrated  that  an  inviscid 
rotational  flow  analysis  of  attached  shock-wave  boundary- layer  interac¬ 
tions  gives  wall  pressure  distributions  which  agree  well  with  experimental 
data  for  turbulent  flows.  In  this  section  calculations  are  presented  using 
the  present  stream  function-vorticity  approach  and  comparisons  are  made 
with  the  asymptotic  theory  results  of  Melnik  and  Grossman  (Ref.  12).  In 
addition  some  observations  are  presented  on  the  occurrence  of  flow 
separation  in  the  application  of  inviscid  rotational  flow  theory.  Finally 
it  IS  demonstrated  that  for  attached  flow,  calculations  using  IBLT  give 
comparable  wall  pressure  distributions  to  that  obtained  from  inviscid 
rotational  flow  theory. 

Inviscid  Rotational  Flow  Analysis 

In  the  application  of  inviscid  rotational  flow  theory  the  upstream 
velocity  profile,  u  (y),  is  chosen  so  as  to  represent  the  outer  portion  of 
a  turbulent  boundary  layer  profile.  The  inner  or  wall  layer  is  ignored.  A 
critical  step  in  the  selection  of  this  profile  is  the  choice  of  magnitude 
of  the  wall  slip  velocity  since  it  is  well  known  that  the  results 
downstream  depend  on  this  value.  In  the  present  calculations,  with  the 
stream  function-vorticity  approach,  it  was  found  that  for  a  given  problem  a 
minimum  value  of  the  slip  velocity  exists.  If  the  slip  velocity  is  speci¬ 
fied  below  this  value  then  inviscid  separation  results  with  the  computed 
wall  pressure  distribution  found  to  depart  significantly  from  that  observed 
experimentally.  This  result  is  easily  explained  by  a  consideration  of  the 
total  and  static  pressure  distributions  as  shown  schematically  in  Fig.  10. 
Here  it  is  seen  that  upstream  of  a  compressive  disturbance,  the  details  of 
which  can  be  ignored,  the  total  pressure  distribution  p^  which  corresponds 
to  the  profile,  u  ,  has  a  minimum  at  the  wall.  This  val&e  of  total 
pressure  remains  constant  along  the  wall  streamline  unless  a  shock  wave 
extends  through  the  rotational  layer  and  impacts  the  surface  and  further 
reduces  the  wall  total  pressure.  If  the  disturbance  is  of  sufficient 
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strength  for  a  given  level  of  inflow  total  pressure,  then  the  wall 
pressure,  as  shown  in  Fig.  10,  will  increase  until  it  is  equal  to  the  wall 
total  pressure  level  thereby  resulting  in  a  stagnation  point  on  the  wall. 
This  point  is  also  a  separation  point  with  the  streamline  pattern  shown  in 
Fig.  10.  Downstream  the  wall  pressure  decreases  due  to  the  finite  reversed 
flow  velocity  and  the  requirement  that  the  total  pressure  remain  constant 
on  the  wall  streamline.  This  expansion  is  probably  not  physically 
meaningful  and  is  a  reflection  of  the  unrealistic  assumption  of  using  an 
inviscid  flow  theory  in  a  separated  flow  region.  A  numerical  example  of 
separated  inviscid  rotational  flow  is  found  in  the  work  of  Taulbee  and 
Robertson*  (Ref.  31).  They  solved  the  incompressible  form  of  the  stream 
function-vorticity  equation  for  the  inviscid  rotational  flow  over  a  forward 
facing  step  on  a  flat  surface.  Their  rotational  flow  model  and  computed 
results  for  two  cases  are  shown  in  Fig.  11.  In  Case  B  the  inviscid  flow 
remains  attached  with  a  monotonic  rise  of  the  wall  pressure  along  the  plate 
ahead  of  the  step;  in  Case  A  inviscid  separation  occurs  near  the  corner 
with  a  dip  in  the  static  pressure  in  that  region. 

In  the  present  rotational  inviscid  flow  analysis  of  transonic  normal 
shock-wave  boundary- layer  interaction,  inviscid  separation  was  avoided  by 
prescribing  the  inflow  wall  slip  velocity  sufficiently  large  so  that  the 
associated  total  pressure  was  greater  than  the  knotm  static  pressure  down¬ 
stream  of  the  shock  wave.  Avoidance  of  inviscid  separation  only  yields  a 
minimum  value  for  the  inflow  slip  velocity;  in  actuality  a  more  specific 
criterion  is  needed  for  subsonic  and  transonic  flows  similar  to  that 
established  by  Rosen,  Roshko,  and  Pavish  (Ref.  32)  for  supersonic  flows 
over  compression  corners. 

With  this  minimum  slip  velocity  criterion  established,  the  discussion 
now  turns  to  a  comparison  of  the  results  obtained  with  the  present  analysis 
with  those  computed  by  Melnik  and  Grossman  (Ref.  12)  from  a  solution  of  the 
transonic  small  disturbance  equation  which  they  derived  from  an  asymptotic 
analysis  of  rotational  flow  interacting  with  a  normal  shock  wave.  Figure 
12  shows  the  upstream  profile  for  the  u-component  of  velocity  which  was 
used  for  the  rotational  inviscid  flow  analysis  of  a  normal  shock-wave, 
boundary-layer  interaction  for  a  free  stream  Mach  number,  M^  1.12.  In 
this  case  the  wall  slip  velocity  was  equal  to  0.68,  which  is  slightly 
higher  than  the  minimum  value  of  0.64  required  to  avoid  flow  separation. 
Figure  12  shows  the  velocity  profile  computed  with  the  present  analysis  at 
the  downstream  boundary,  x  ~  15;  it  is  observed  that  a  large  decrease 
occurs  in  u  in  the  wall  region  which  would  lead  to  separation  if  the  inflow 
slip  velocity  was  specified  below  the  minimum  value.  The  upstream  profile 
used  by  Melnik  and  Grossman,  which  is  also  shown  in  Fig.  12,  extends  to  the 
wall  without  regard  to  the  minimum  slip  velocity  criterion  discussed  above. 
In  their  small  disturbance  theory  the  static  pressure  rises  monotonically 
to  the  imposed  downstream  level  independent  of  the  occurrence  of  flow 
reversal  and  total  pressure  considerations  since  the  pressure  increase  is 
assumed  to  be  proportional  to  the  velocity  decrease,  that  is,  Ap  ~  -Au. 

The  present  computations  were  made  in  the  region  -13  <  x  <  15  and  0  < 
y  <  8  with  the  normal  shock  imposed  at  x  *  0.  Recall  that  all  lengths  are 
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nondinensionalized  by  the  upstream  boundary-layer  thickness.  In  order  to 
facilitate  the  comparison  of  the  present  results  with  those  of  Melnik  and 
Grossman  the  approximate  shock  jump  condition  which  they  imposed  was  used 
in  the  present  calculation.  This  condition  gives  a  velocity  decrease  of  Lu 
■  -0.212  and  a  corresponding  pressure  increase  &p  ~  0.212,  in  contrast  to 
the  exact  Rankine-Hugoniot  values  of  Au  ■  -0.169  and  Ap  “  0.169  for  * 
1.12.  The  present  results  were  obtained  with  the  zebroid  line  relaxation 
scheme  which  was_|^epeated  until  the  maximum  difference  equation  residual 
was  less  than  10  .  The  vorticity  and  entropy  were  updated  after  every  10 

sweeps  of  the  flow  field.  In  these  calculations  81  grid  points  were  evenly 
distributed  in  the  x-direction  which  resulted  in  a  streamwise  step  size  of 
Ax  >  0.375.  In  the  y-direction  21  grid  points  were  distributed 
nonuniformly  with  a  coordinate  stretching  used  to  concentrate  10  of  these 
grid  points  in  the  boundary  layer,  0  <  y  <  1.  These  grid  distributions 
differ  from  that  used  by  Melnik  and  Grossman  as  they  solved  in  the  infinite 
half-plane,  -“<x<“^,  0<y<°=  with  coordinate  stretchings  used  in  both 
the  X  and  y  directions  to  concentrate  the  mesh  (78  by  577  in  the  boundary 
layer  and  near  the  shock  wave.  Melnik  and  Grossman  reported  that  various 
studies  were  made  on  the  effect  of  the  mesh  on  the  accuracy  of  the  method; 
a  similar  study  should  be  made  in  the  future  with  the  present  calculations 
to  determine  the  sensitivity  of  the  results  to  the  mesh  size  and  distribu¬ 
tion.  One  test  that  was  performed  in  the  present  study  was  to  move  the 
the  outer  computational  boundary  from  y  »  8  to  y  *  16.  Only  a  small 
difference  was  observed  in  the  computed  wall  pressure  distribution. 

Figure  13  shows  a  comparison  of  the  wall  pressure  obtained  from  the 
present  formulation  of  inviscid  rotational  flow  in  comparison  with  that 
predicted  by  the  asymptotic  small  disturbance  theory  of  Melnik  and 
Grossman.  The  good  agreement  observed  between  these  results  is  encouraging 
as  there  are  substantial  differences  in  the  two  approaches.  Later  after 
Melnik  and  Grossman  first  presented  this  result,  they  included  the  effect 
of  boundary  layer  blockage  in  an  axisymmetric  channel  so  that  they  could 
compare  their  predicted  wall  pressure  with  the  experimental  pressure 
distribution  measured  by  Gadd  (Ref.  17)  for  this  M^  >  1.12  case.  They 
obtained  good  agreement  with  the  data  thereby  verifying  the  accuracy  of 
inviscid  rotational  flow  theory  for  attached  transonic  normal  shock-wave, 
boundary- layer  interaction. 

Before  ending  this  section,  it  is  pointed  out  that  it  was  found 
imperative  in  the  present  stream  function-vorticity  method  to  precisely 
satisfy  the  Crocco  relation,  Eq.  7.  Initially,  calculations  were  made 
using  the  parallel  flow  approximation  in  which  the  vorticity  to  pressure 
ratio  and  the  entropy  were  assumed  constant  along  lines  of  y  *  constant 
instead  of  along  streamlines.  Although  in  this  case  the  streamlines  only 
depart  slightly  from  constant  ylines,  nonetheless  a  slight  error  in  the 
entropy  a  is  exponentially  amplified  in  the  computation  of  the  pressure  in 
Eq.  6  which  led  to  spurious  pressure  variations  across  the  boundary  layer 
in  the  present  calculations.  These  problems  were  eliminated  by  exactly 
satisfying  the  Crocco  equation  with  a  curve  fit  of  the  entropy-stream 
function  distribution  at  the  inital  station  and  the  use  of  analytical 
differentiation  to  obtain  the  vorticity  to  pressure  ratio  at  any  desired 
point  in  the  flow  field. 
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Interacting  Boundary  Laver  Theory  (IBLT) 

The  discussion  in  the  previous  section  demonstrates  that  a  continuous 
pressure  rise  is  produced  at  the  wall  if  the  outer  layer  of  a  turbulent 
boundary  layer  is  modeled  as  an  inviscid  rotational  layer  which  interacts 
with  a  normal  shock  wave.  A  question  which  has  been  briefly  addressed  in 
the  present  study  is  whether  or  not  IBLT  is  capable  of  capturing  the  proper 
interaction  scales  of  an  attached  transonic  shock-wave,  boundary- layer 
interaction.  Only  a  first-order  analysis  has  been  made  for  this  case  of 
Moo  ■  1.12  since  it  was  found,  as  will  be  discussed  later,  that  inclusion  of 
normal  pressure  gradients  in  the  viscous  layer  for  a  much  stronger  shock 
wave  with  >>  1.4  had  only  a  minor  impact  on  the  wall  pressure,  displace¬ 
ment  thickness,  and  skin  friction  distributions. 

In  the  IBLT  calculation,  the  inviscid  computation  was  made  on  the  same 
grid  discussed  previously.  Uniform  inviscid  inflow  was  now  assumed  at  the 
upstream  boundary  which  results  in  irrotational  inviscid  flow  throughout 
the  computational  region  since  the  entropy  rise  at  the  shock  wave  was 
neglected  in  the  present  investigation.  The  viscous  analysis,  which  is  a 
first-order  inverse  boundary- layer  calculation,  was  computed  on  the  same 
streamwise  mesh  as  the  inviscid  analysis  with  100  grid  points  distributed 
nonuniformly  across  the  boundary  layer.  The  upstream  profile  was  the  same 
as  that  used  in  the  inviscid  rotational  flow  calculation  only  now  modified 
near  the  wall  to  satisfy  the  no  slip  condition.  The  IBLT  calculation  was 
made  with  a  relaxation  factor  of  u  ■■  .75  in  the  update  formula  in  £q.  (68); 
30  global  iterations  were  required  to  reduce  the  difference  betwee^  the 
viscous  and  inviscid  predictions  for  the  edge  velocity  down  to  10 

Figure  13  shows  the  computed  wall  pressure  from  IBLT  in  comparison 
with  the  inviscid  rotational  flow  results  obtained  by  Melnik  and  Grossman 
and  the  present  theory.  The  smooth  pressure  distribution  obtained  with 
IBLT  demonstrates  that  displacement  thickness  interaction  results  in  an 
interaction  length  consistent  with  that  obtained  from  inviscid  rotational 
flow  theory  for  a  transonic  normal  shock  wave  impinging  on  a  turbulent 
boundary  layer.  Figure  14  shows  the  computed  displacement  thickness  and 
skin  friction  distribution  obtained  for  this  case  from  IBLT.  As  expected 
the  flow  remains  attached  in  this  case  as  separation  generally  is  thought 
to  occur  for  Mach  numbers  greater  than  about  1.3.  Nonetheless  the  present 
results  show  that  even  in  this  case  the  boundary- layer  response  is  not 
trivial  as  the  skin  friction  coefficient  is  reduced  in  half  and  the 
displacement  thickness  doubles  in  the  small  streamwise  length  of  this 
interaction. 

It  is  observed  in  Fig.  13  that  the  streamwise  length  scale  of  the 
overall  wall  pressure  rise  is  only  about  10  upstream  boundary- layer 
thicknesses.  Since  the  shock  wave  in  the  inviscid  flow  is  computed  with  a 
shock  capturing  technique,  there  was  initial  concern  that  the  thickness  of 
the  shock  transition  region  induced  by  the  numerical  technique  could  be 
comparable  to  that  caused  by  the  actual  boundary  layer.  Figure  15  shows 
the  wall  pressure  obtained  from  a  numerical  experiment  which  was  performed 
to  assess  the  numerical  smearing  of  the  shock  wave.  In  this  calculation 
the  effect  of  the  boundary  layer  was  removed  and  only  a  uniform,  irrota¬ 
tional  inviscid  flow  subject  to  a  discontinuity  at  the  outer  computational 
boundary  was  computed.  The  wall  pressure  in  Fig.  15  shows  that  although 
the  shock  is  slightly  smeared  at  the  wall  the  streamwise  length  scale  of 


this  pressure  rise  is  significantly  less  than  that  induced  by  the  actual 
boundary  layer.  Similar  tests  to  this  were  performed  throughout  the 
present  investigation  which  verified  that  the  artificial  viscosity  used  in 
the  inviscid  scheme  had  a  minimum  impact  on  the  overall  results. 

Separated  Flow  Calculations 

Attention  is  now  focused  on  the  analysis  of  transonic  normal  shock- 
wave,  turbulent  boundary-layer  interactions  in  which  the  shock  wave  is  of 
sufficient  strength  to  separate  the  boundary  layer.  The  detailed  experi¬ 
mental  data  of  Rooi  (Refs.  5  and  6)  has  been  extensively  used  in  this 
investigation  to  evaluate  the  present  IBLT  calculations  with  and  without 
the  inclusion  of  normal  pressure  gradients.  In  this  section  both  free  air 
and  bounded  channel  calculations  are  presented  for  the  Rooi  data.  Also, 
inviscid  rotational  flow  calculations  are  compared  with  Rooi's  wall 
pressure  distribution  for  a  case  in  which  flow  separation  occurs. 

Free  Air  Results 

Rooi  performed  a  detailed  experimental  investigation  of  transonic 
shock-wave,  boundary- layer  interaction  in  the  channel  shown  schematically 
in  Fig.  16.  A  normal  shock  wave  was  established  in  this  passage  by  placing 
a  "choking  wedge"  on  the  lower  flat  plate;  the  choking  wedge  was  positioned 
in  the  stream  direction  such  that  a  normal  shock  wave  emanated  from  the 
leading  edge  of  the  upper  flat  plate.  This  shock  wave  impinged  on  the 
lower  flat  plate  boundary  layer  which  was  tripped  to  turbulent  near  the 
leading  edge  as  shown  in  Fig.  16.  For  *  1.4,  which  is  the  case  computed 
in  the  present  investigation,  Rooi  experimentally  deduced  that  the 
boundary- layer  thickness  6*  on  the  lower  surface,  just  upstream  of 
interaction  region,  is  6.6  mm.  This  value  of  6^'  has  been  used  to  nondimen- 
sionalize  all  lengths  in  the  present  results  and  in  the  presentation  of  the 
experimental  measurements  of  Rooi.  We  have  used  Rooi’s  data  exclusively 
from  Ref.  5;  later  he  investigated  the  effects  of  higher  Mach  numbers  by 
testing  at  >  1.44  and  1.46. 

An  attempt  was  made  to  compute  the  Rooi  >  1.4  case  with  the  outer 
computational  boundary  treated  as  a  "free  air"  boundary  as  shown  in  Fig.  3. 
This  "free  air"  boundary  was  placed  at  the  location  of  the  upper  wall  shown 
in  Fig.  16  which  is  at  y  >  12.1  (80  mm).  Along  this  boundary  ip  was  set 
equal  to  a  constant  with  p  and  u  set  equal  to  1  ahead  of  the  shock  wave  and 
at  the  respective  Rankine-Hugoniot  values  downstream  of  the  shock.  The 
z-coordinate  has  been  redefined  as  z  ■  (x  -z*)/6*  where  x*  is  the 
z-location  of  the  start  of  interaction.  A  vSlue^of  x*  >  ^26  mm  has  been 
inferred  from  Rooi's  measurements  (Ref.  5)  and  has  been  used  to  scale  the 
z-coordinate  both  in  the  theoretical  and  experimental  results. 

The  free  air  IBLT  calculation  of  the  Rooi  -  1 .4  case  was  made  in 
the  computational  region,  -4  <  z  <  26  and  0  <  y  <  12.  In  the  inviscid 
calculation  100  grid  points  were  evenly  spaced  in  the  z-direction  with  Ax  * 
.30;  in  the  y-direction  a  coordinate  transformation  was  used  to  establish  a 
nonuniform  grid  distribution  with  the  smallest  spacing  of  ^y  «  .1025  placed 
at  the  wall.  In  the  viscous  calculation  the  streamwise  grid  was  the  same 
as  that  used  in  the  inviscid  calculation;  in  the  y-direction  100  grid 
points  were  nonuniformly  distributed  across  the  boundary  layer  with  the 


smallest  spacing  placed  near  the  vail.  The  viscous  grid  was  extended  with 
a  constant  9  percent  stretching  from  the  wall  to  a  maximum  n  of  22.  In 
this  calculation  the  effects  of  normal  pressure  gradients  have  not  been 
included  and  hence  n  *  n  with  h  -  1.  A  total  of  40  global  viscous-inviscid 
iterations  were  made  to  obtain  a  converged  solution. 

Figure  17  shows  the  present  IBLT  calculation  for  the  wall  pressure, 
displacement  thickness,  and  skin  friction  in  comparison  with  the 
experimental  measurements  of  Kooi  (Ref.  5).  The  Reynolds  number  for  this 
flow  is  2.97  X  10^J>ased  on  the  free  stream  quantities  and  the  boundary- 
layer  thickness,  6  -  6.6  mm,  upstream  of  the  interaction.  Figure  17 

shows  significant  differences  between  the  theoretical  and  experimental 
results;  these  differences  are  primarily  due  to  the  boundary- layer  blockage 
effects  which  existed  in  the  experiment  that  prevented  the  experimental 
wall  pressure  from  rising  downstream  to  the  Rankine-Hugoniot  normal  shock 
value.  In  the  experiment  the  large  boundary- layer  growth  produced  by  the 
shock  wave  accelerates  the  subsonic  channel  flow  downstream  of  the  shock 
wave  and  therefore  has  a  significant  effect  on  the  interaction  flow  field 
results  as  shown  in  Fig.  17.  In  contrast  the  IBLT  calculation  shows  that, 
if  the  pressure  rises  to  the  Rankine-Hugoniot  level,  then  the  overall 
growth  of  the  displacement  thickness  and  the  streamwise  extent  of  the 
separated  flow  region  is  much  larger  than  that  measured  experimentally. 
Fortunately  we  were  able  to  account  for  the  blockage  effect  in  this  experi¬ 
ment  by  imposing  in  the  IBLT  calculation  the  static  pressure  distribution 
measured  by  Kooi  along  the  line  y  «  45  mm  or  y  >  6 .8  above  the  lower  flat 
plate.  This  computational  region  located  in  the  interior  of  the  Kooi 
channel  is  shorn  in  Fig.  16  and  will  be  utilized  in  the  calculations  to  be 
discussed  in  the  remainder  of  this  report.  Before  the  IBLT  results  are 
presented  a  brief  examination  is  made  of  the  use  of  inviscid  rotational 
flow  theory  for  the  prediction  of  separated  transonic  shock-wave, 
boundary- layer  interaction. 

Inviscid  Rotational  Flow  Analysis 

Figure  18  shows  the  inflow  profiles  of  velocity  and  the  corresponding 
entropy-stream  function  distribution  used  in  the  inviscid  rotational  flow 
analysis  of  the  Kooi  ■  1.4  case.  The  velocity  profile  was  prescribed  so 
that  it  matches  the  experimental  data  over  most  of  the  boundary  layer; 
however  due  to  the  minimum  slip  velocity  constraint,  which  was  previously 
discussed,  the  prescribed  inflow  velocity  deviates  from  the  experiment  near 
the  wall  as  shown  in  Fig.  18.  Also  shown  in  Fig.  18  is  a  prediction  of 
this  initial  boundary  layer  profile  using  the  first-order  direct  version  of 
the  present  boundary- layer  analysis  subject  to  a  zero  pressure  gradient. 
This  calculation  was  started  at  the  leading  edge  of  the  flat  plate  with 
instantaneous  transition  from  laminar  to  turbulent  flow  assumed  to  occur  at 
the  location  of  the  experimental  trip  strip,  shown  in  Fig.  16.  Excellent 
agreement  is  shown  here  between  the  predicted  and  measured  results. 

Figure  19  shows  the  computational  mesh  used  in  this  inviscid 
rotational  flow  analysis.  The  mesh  consists  of  81  evenly  distributed 
x-grid  points  with  Ax  *  .30  and  21  y-grid  points  distributed  as  shown  in 
Fig.  19.  Approximately  8  points  were  placed  in  the  boundary  layer  with  the 
minimum  Ay  >  .0883  located  at  the  wall. 
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A  comparison  of  the  wall  pressure  computed  from  inviscid  rotational 
flow  theory  with  Kooi's  experimental  measurements  is  shown  in  Fig.  20.  In 
contrast  with  the  inviscid  rotational  theory  result  shown  in  Fig.  13  for 
>  1.12,  the  present  calculation  at  =  1.4  shows  a  nearly  discontinuous 
rise  in  the  wall  pressure  at  the  imposed  shock  location  of  x  <=3.6.  The 
lack  of  upstream  influence  is  expected  in  the  =<  1.4  case  since  the 
entire  upstream  inviscid  flow  was  supersonic  with  the  minimum  Mach  number, 

M  >1.13,  occurring  at  the  wall.  In  the  >  1.12  case  the  sonic  line  of 
tKe  entering  profile  was  located  at  y  >  .4  with  a  subsonic  slip  Mach  number 
of  M  >  .72.  It  should  be  noted  that  the  lack  of  upstream  influence  does 
not  prevent  inviscid  rotational  flow  theory  from  being  successfully  used  to 
describe  supersonic  interacting  flow  fields.  For  example,  Roshko  and 
Thompke  (Ref.  33)  showed  very  good  agreement  between  theory  and  experiment 
for  the  M^  <=  3.93  flow  over  a  15  deg  compression  corner.  In  this  case  the 
flow  remained  attached  with  a  discontinuous  rise  in  the  wall  pressure  at 
the  corner  followed  downstream  by  a  continuous  compression  to  the  wedge 
pressure  level. 

The  poor  agreement  shown  in  Fig.  20  demonstrates  the  inadequacy  of 
inviscid  rotational  flow  theory  for  the  prediction  of  the  wall  pressure 
distribution  when  the  shock  wave  is  of  sufficient  strength  to  separate  the 
boundary  layer.  In  contrast  with  the  attached  flow  case  where  the  viscous 
wall  layer  is  passive  and  hence  can  be  ignored  the  increase  in  shock 
strength  now  results  in  such  a  large  increase  in  pressure  that  the  wall 
layer  grows  significantly  thereby  displacing  the  outer  portion  of  the 
boundary  layer.  Therefore  the  wall  layer  must  be  actively  included  in  the 
interaction  calculation.  In  the  two-layer  IBLT  technique,  for  which 
separated  flow  results  are  presented  in  the  next  section,  the  viscous  layer 
is  treated  in  its  entirety  without  a  breakdown  into  an  outer  inviscid 
rotational  layer  and  a  viscous  wall  layer.  This  latter  viscous-inviscid 
interaction  flow  model  has  been  designated  as  a  three-layer  model  as  shown 
in  Fig.  2  and  is  presented  in  Appendix  A. 

Channel  Results  for  IBLT  with  Normal  Pressure  Gradients 

In  this  section  comparisons  are  presented  between  the  present 
computations  of  IBLT  with  normal  pressure  gradients  and  the  transonic  shock 
induced  separation  data  of  Kooi  (Ref.  5).  First  the  details  of  the 
numerical  calculation  will  be  given  and  then  a  brief  discussion  of  a 
modification  to  the  Cebeci-Smith  algebraic  turbulence  model  which  was  found 
to  improve  the  agreement  of  the  theoretical  prediction  with  the  data. 

The  IBLT  calculations  were  performed  in  the  "channel"  computational 
region  shown  in  Fig.  16  with  the  experimental  pressure  distribution  (shown 
in  Fig.  23)  imposed  along  the  outer  boundary  (y  >  6.8)  to  account  for  the 
boundary- layer  blockage  effect  which  was  previously  discussed.  The 
inviscid  grid  distribution  that  was  used  is  the  same  as  that  shown  in  Fig. 
19  which  was  used  in  the  rotational  inviscid  calculation.  The  inclusion  of 
normal  pressure  gradients  in  the  IBLT  calculation  requires  a  fine  grid  to 
be  used  in  the  inviscid  analysis  near  the  surface  to  provide  an  accurate 
description  of  the  equivalent  inviscid  flow  (EIF)  quantities  which  are 
input  to  the  viscous  layer  analysis.  The  viscous  calculation  used  the  same 
streamwiae  mesh  as  the  inviscid  analysis  to  avoid  interpolation  which  has 
been  found  in  previous  interaction  calculations  (Ref.  30)  to  cause  oscilla¬ 
tions.  In  the  normal  direction,  101  grid  points  were  nonuniformly 
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distributed  from  the  wall  across  the  viscous  layer  to  n  *  1.45.  The 
iteration  between  the  viscous  and  inviscid  analyses  were  performed  with  tne 
global  technique  shown  in  Fig.  9.  In  the  present  calculations  a  converged 
solution  was  obtained  first  for  first-order  IBLT  in  which  the  pressure  was 
assumed  constant  across  the  boundary  layer.  This  calculation  in  which  Z  = 
.4,  required  SO  global  iterations  to  reduce  t^e  difference  between  the 
viscous  and  inviscid  prediction  for  u^  to  10  .  These  results  were  then 

used  to  initialize  the  IBLT  calculation  with  normal  pressure  gradients.  In 
the  latter  calculation  it  was  found  necessary  to  reduce  ti  to  .2  to  obtain 
convergence  which  was  achieved  in  typically  20  global  iterations.  The 
reduction  of  the  global  relaxation  factor  is  not  surprising  since  the 
coupling  between  the  viscous  and  inviscid  flows  is  significantly  greater 
with  the  inclusion  of  the  numerical  results  for  the  equivalent  inviscid 
flow  in  the  viscous  defect  equations. 

It  is  well  known  that  use  of  the  Cebeci-Smith  algebraic  turbulence 
model  generally  results  in  an  overprediction  of  the  pressure  and  an  under- 
prediction  of  the  displacement  thickness  in  the  analysis  of  separated 
turbulent  flow.  In  the  past  it  has  been  found  that  some  improvement  in  the 
predicted  results  could  be  obtained  by  a  reduction  in  the  local  values  of 
the  eddy  viscosity  coefficient  predicted  by  this  model.  One  example  of 
this  type  of  modification  is  the  Shang-Hankey  (Ref.  34)  relaxation  model 
which  has  been  used  by  a  number  of  investigators  for  the  computation  of 
separated  turbulent  flow.  In  this  model  both  the  inner  and  outer  layer 
eddy  viscosity  levels  are  constrained  to  depart  by  only  a  specified  amount 
from  an  upstream  equilibrium  boundary  layer  level.  More  recently  a  simple 
numerical  experiment  by  Carter  (Ref.  35)  demonstrated  that  a  reduction  in 
the  outer  region  length  scale,  which  Clauser  (Ref.  36)  deduced  as  the 
incompressible  displacement  thickness  for  mild  pressure  gradient  flows, 
resulted  in  similar  improvements  to  the  separated  flow  predictions  as 
demonstrated  by  others  with  the  Shang-Hankey  model.  In  a  recent  study  con¬ 
ducted  by  Edwards  (Ref.  37),  the  incompressible  turbulent  separated  flow 
data  of  Simpson,  et  al.  (Ref.  38)  was  examined  to  determine  if  a  modified 
outer  length  scale  could  be  inferred  from  their  detailed  measurements  for 
the  streamwise  velocity  profiles  and  the  Reynolds  stress  distributions.  It 
was  found  that  if  the  outer  eddy  viscosity  coefficient  was  based  on  the 
incompressible  momentum  thickness  instead  of  the  incompressible 
displacement  thickness  then  a  significant  improvement  was  found  in  the 
correlation  of  the  maximum  Reynolds  shear  stress  with  the  local  gradient  of 
the  mean  velocity  profile.  With  this  modification  the  outer  region  eddy 
viscosity  coefficient  is  expressed  as 

y  =0.0168  Hp  uee,nc  (69) 

where  the  quantity  H  is  the  flat  plate  incompressible  shape  factor 
6*  /6.  equal  to  9/7,  and  the  quantity  6.  is  the  incompressible 

mraentm^ thickness.  For  flat  plate  flows  ISis  model  is  identical  to 
Clauser 's  model.  As  will  be  seen  in  this  section  this  modified  Cebeci- 
Smith  (CS)  turbulence  model  has  improved  the  present  predictions  for 
transonic  shock  wave,  boundary- layer  interaction. 

Figure  21  shows  a  comparison  of  the  present  computations  from  IBLT 
with  the  experimental  data  of  Kooi  (Ref.  5)  for  the  wall  pressure, 


displacement  thickness,  skin  friction  coefficient,  momentum  thickness  and 
ahape  tactor  for  >  1.4.  Kooi  used  the  same  generalized  definitions  of 
the  displacement  thickness  and  the  momentum  thickness  in  presenting  his 
experimental  results  as  has  been  used  in  the  present  theoretical  study. 

The  skin  friction  coefficient  is  based  on  the  flow  properties  at  the 
boundary- layer  edge  which  in  the  calculations  has  been  defined  where  F  = 
.995F.  and  in  the  experiment  as  the  point  where  the  pitot  pressure  became 
constant.  The  computed  results  are  shown  for  both  the  CS  and  the  modified 
CS  turbulence  model  discussed  above.  The  position  of  the  imposed  normal 
shock  wave  is  indicated  in  Fig.  21.  Figures  21(a)  and  (b)  show  that  good 
agreement  is  obtained  with  the  data  with  the  present  IBLT  using  the  modi¬ 
fied  CS  turbulence  model.  The  agreement  shown  for  the  wall  pressure  is 
substantially  better  than  that  shown  earlier  for  the  free  air  IBLT  calcu¬ 
lation  and  that  predicted  by  the  inviscid  rotational  flow  analysis.  Com¬ 
parison  of  Figs.  21(a)  and  (b)  show  that  the  turbulence  model  modification 
simultaneously  improves  the  agreement  for  the  wall  pressure  and  the 
displacement  thickness  which  substantiates  the  use  of  the  displacement  body 
concept  in  separated  flows.  The  distributions  of  the  skin  friction 
coefficient  in  Fig.  21(c)  compare  the  extent  of  the  separated  region  in  the 
present  calculations  with  that  deduced  by  Kooi  from  his  measured  velocity 
profiles.  Both  of  the  results  obtained  with  the  CS  and  modified  CS 
turbulence  model  show  good  agreement  with  the  data  in  the  approach  to  the 
separation  point  which  occurs  just  upstream  of  the  location  of  the  imposed 
normal  shock  wave.  However,  downstream  the  modified  model  underpredicts 
the  experimental  skin  friction  and  therefore  predicts  a  much  larger 
separated  flow  region  than  that  found  experimentally  or  in  the  present 
calculation  using  the  unmodified  CS  turbulence  model.  Modification  to 
the  outer  region  of  the  CS  turbulence  model  has  resulted  in  an  improved 
prediction  of  the  wall  pressure  and  displacement  thickness;  the  lack  of 
agreement  for  the  skin  friction  indicates  that  further  work  is  required  to 
correspondingly  improve  the  accuracy  of  the  IBLT  with  the  modified  CS 
turbulence  model  in  the  near  wall  region.  A  similar  trend  is  observed  in 
Fig.  21(d)  for  the  momentum  thickness  comparisons  in  which  it  is  seen  that 
the  aiodified  model  underpredicts  the  measured  6  distribution.  This  result 
is  not  surprising  since  the  momentum  thickness,  as  well  as  the  skin 
friction,  is  strongly  influenced  by  the  near  wall  solution.  Finally,  the 
distributions  of  the  shape  factor,  H  •=  6*/e,  are  shown  in  Fig.  21(e)  in 
which  it  is  seen  that  neither  model  accurately  predicts  the  experimental 
distribution. 

Figure  22  shows  a  comparison  of  the  computed  and  the  experimental 
profiles  for  the  u-component  of  velocity  at  7  streamwise  locations  upstream 
and  downstream  of  the  normal  shock  wave.  The  insert  in  Fig.  22  shows  the 
X- location  of  these  profiles;  the  corresponding  station  number  used  by 
Kooi  (Ref.  5)  to  label  these  profiles  has  also  been  included.  Linear 
streasiwise  interpolation  was  used  in  the  calculations  to  obtain  the 
computed  profiles  at  the  same  x-location  as  that  of  the  experimental 
measurements.  The  edge  velocity,  u  ,  used  to  normalize  the  u-profile  was 
determined  by  the  same  criterion  mentioned  earlier  for  the  respective  theo¬ 
retical  and  experimental  boundary- layer  thicknesses. 
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The  velocity  profiles  shown  in  Fig.  22  illustrate  the  rapid  changes 
which  occur  in  a  turbulent  boundary- layer  profile  due  to  an  incident  normal 
shock  wave.  For  example,  comparison  of  the  experimental  profile  at  x  ^  .4 
with  that  just  downstream  of  the  shock  wave  at  x  =3.8  shows  the  enormous 
change  which  occurs  in  only  3  upstream  boundary- layer  thicknesses.  The 
velocity  profiles  in  Fig.  22  show  that  the  IBLT  results  obtained  with  the 
modified  CS  turbulence  model  are  in  overall  better  agreement  with  the 
experimental  data  than  those  results  obtained  without  this  modification. 

The  agreement  is  particularly  good  in  the  outer  region  of  the  turbulent 
boundary  layer  where  this  modification  was  made.  Downstream  of  the  shock 
wave  the  modified  model  underpredicts  the  velocity  recovery  in  the  near 
wall  region  as  was  expected  from  the  skin  friction  and  momentum  thickness 
comparisons  shown  in  Fig.  21.  The  results  obtained  with  the  unmodifie 
model  show  a  higher  velocity  in  the  wall  region  at  each  location  than  the 
modified  model  which  is  consistent  with  the  skin  friction  difference  shown 
in  Fig.  21.  Despite  the  complexity  of  this  flow,  it  is  observed  in  Fig.  22 
that  the  overall  growth  of  the  viscous  layer,  which  appears  to  have 
distinct  edge,  is  well  predicted  by  the  present  IBLT  with  the  modified  CS 
turbulence  model.  The  details  of  the  flow  at  the  edge  of  the  viscous  layer 

where  the  real  viscous  flow  (RVF)  and  the  equivalent  inviscid  flow  (EIF) 

merge  will  be  discussed  later. 

Figure  23  shows  a  comparison  of  the  computed  and  measured  static 

pressure  distributions  at  various  y- levels  in  the  flow  field.  The  IBLT 

results  which  are  shown  are  those  obtained  with  the  modified  CS  turbulence 
model.  As  mentioned  previously,  the  static  pressure  distribution  at  y  = 
6.8,  which  is  the  location  of  the  outer  inviscid  computational  boundary, 
was  set  equal  to  that  measured  by  Kooi  in  order  to  account  for  the  effect 
of  boundary  layer  blockage.  At  the  other  y- locations  shown  in  Fig.  23  the 
agreement  of  the  IBLT  result  with  the  experiment  is  good  except  just 
downstream  of  the  shock  wave  where  the  pressure  is  overpredicted.  Further 
work  is  needed  to  eliminate  this  pressure  overshoot.  Figure  23  shows  that 
the  static  pressure  varies  between  the  wall  and  y  =  2.3.  Later  in  this 
section  the  detailed  pressure  variation  across  the  viscous  layer  will  be 
presented. 

Figure  24  shows  the  flow  structure  deduced  from  the  present  calcula¬ 
tion  with  the  modified  CS  turbulence  model  in  comparison  with  that  deduced 
by  Kooi  from  his  measurements  and  Schlieren  photographs.  Outside  of  the 
immediate  shock  wave  regions,  the  theoretical  and  experimental  Mach 
contours  in  the  inviscid  flow  region  show  the  same  pattern,  particularly 
do%rastream  of  the  shock  wave.  Overall  the  agreement  between  the  predicted 
and  measured  flow  structure  is  considered  reasonable.  Further  research  is 
needed  to  determine  why  the  triple  or  lambda  shock  structure  observed  in 
the  experiment  is  not  present  in  the  theoretical  flow  pattern.  It  is 
possible  that  this  difference  is  due  to  the  neglect  of  the  entropy  rise  or 
total  pressure  decrease  across  the  shock  in  the  IBLT  calculation. 

The  variation  of  the  pressure  across  the  boundary  layer  from  the  equi¬ 
valent  inviscid  flow  (EIF)  is  sho%m  in  Fig.  25  for  selected  streamwise  sta¬ 
tions  both  upstream  and  do%rastream  of  the  shock  wave  which  is  located  at 
X  =  3.6.  The  x-location  of  each  of  these  normalized  pressure  profiles  is 
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shown  in  a  scaled  drawing  inserted  into  Fig*  25.  Since  the  variation  in 
the  static  pressure  across  the  viscous  layer  was  insensitive  to  the  turbu¬ 
lence  model  only  the  IBLT  results  with  the  modified  CS  turbulence  model 
are  shown  in  Fig.  25.  The  arrows  on  each  of  the  profiles  are  the  edge  of 
the  viscous  layer.  The  kink  in  the  profiles  rear  the  wall  is  at  the  dis¬ 
placement  thickness  since,  as  was  discussed  earlier,  the  pressure  was 
assumed  constant  from  6*  to  the  wall.  Figure  25  shows  the  anticipated 
trend  in  the  static  pressure  variation  in  that  in  the  immediate  vicinity  of 
the  shock  wave  the  pressure  decreases  with  increasing  y  upstream  of  the 
shock  wave  and  increases  with  increasing  y  downstream  of  the  shock  wave. 

The  variation  of  the  static  pressure  across  the  viscous  layer  reaches  a 
maximum  of  about  10  percent  in  the  immediate  vicinity  of  the  impinging 
shock  wave;  as  expected,  this  variation  vanishes  upstream  and  downstream  of 
this  region  as  the  flow  returns  to  that  of  a  conventional  constant  pressure 
boundary  layer . 


Figures  26  shows  the  u-velocity  component  functions  from  the  real 
viscous  flow  (RVF)  and  equivalent  inviscid  flow  (EIF)  at  the  same 
streamwise  stations  both  upstream  and  downstream  of  the  shock  wave  as  were 
shown  for  the  viscous  layer  pressure  in  Fig.  25.  The  EIF  and  RVF  stream 
functions  are  eho%m  in  Fig.  27  at  several  streamwise  locations.  As 
discussed  previously  for  the  static  pressure,  tj^e  EIF  for  the  u-velocity 
and  density  were  assumed  constant  for  0  <  y  <  6  in  the  present 


calculations.  Hence  the  stream  function  varies  linearly  in  this  region. 
Figures  26  and  27  show  the  smooth  merging  between  the  RVF  and  EIF  as  the 
edge  of  the  viscous  layer  is  approached.  No  numerical  difficulties  were 
encountered  in  the  implementation  of  these  outer  boundary  conditions. 
Examination  of  the  velocity  profiles  show  the  same  trends  as  that  observed 
for  the  static  pressure,  that  is,  upstream  of  the  shock  wave  the  decrease 
in  the  pressure  from  the  wall  to  the  edge  the  viscous  layer  results  in  a 
correspnding  velocity  increase  as  shown  for  the  EIF  at  x  2.0,  2.9  and 

3.5.  Downstream  of  the  shock  wave  the  opposite  trend  occurs  as  evidenced 
by  the  decrease  in  velocity  across  the  viscous  layer  in  the  EIF  for  x  > 

3.6.  Then  further  downstream  as  seen  in  Fig.  26(c)  the  velocity  from  the 
EIF  varies  much  less  across  the  boundary  layer  as  the  strong  interaction 
between  the  viscous  and  inviscid  flows  subsides.  These  results  demonstrate 
the  present  capability  to  account  for  the  variation  of  the  inviscid  flow  on 
the  scale  of  the  boundary  layer;  however  the  variation  of  the  EIF  in  the 
viscous  layer  is  small  as  seen  in  the  small  changes  in  the  pressure, 
velocity  and  the  linearity  of  the  stream  function  in  this  region.  This 
suggests  that  the  inclusion  of  normal  pressure  gradients  play  a  secondary 
role  in  the  present  calculations.  This  topic  is  addressed  in  further 
detail  in  the  last  section  of  this  report. 


The  viscous  streamline  pattern  computed  for  this  «  1.4  case  is 
shown  in  Fig.  28.  The  vertical  scale  has  been  enlarged  in  order  to  display 
the  details  of  the  separation  bubble.  Only  the  streamlines  below  the  dis¬ 
placement  thickness  are  shown.  It  is  observed  in  Fig.  28  that  separation 
occurs  just  upstream  of  the  normal  shock  wave  with  a  rapid  turn  in  the 
viscous  streamlines  propagating  across  the  viscous  layer.  Downstream  of 
the  Shock  wave  the  streamline  pattern  is  smoother.  The  height  of  the 
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separation  bubble  increases  to  about  0.1  of  the  upstream  undisturbed 
boundary-layer  thickness. 


Comparison  of  IBLT  With  and  Without  Normal  Pressure  Gradient 

One  of  the  advantages  of  the  present  approach  is  that  IBLT 
computations  can  be  made  with  and  without  normal  pressure  gradients  in 
order  to  assess  the  relative  importance  of  including  pressure  variations 
across  the  viscous  layer.  With  the  pressure  assumed  constant  across  the 
viscous  layer  the  present  analysis  reverts  back  to  a  first-order  inverse 
boundary  layer  scheme  interacting  with  the  outer  inviscid  flow  through  dis¬ 
placement  thickness  coupling.  Figure  29  shows  a  comparison  of  the  wall 
pressure,  displacement  thickness,  and  skin  friction  computed  from  IBLT  with 
and  without  normal  pressure  gradients.  These  calculations  were  made  with 
the  modified  CS  turbulence  model.  Comparison  of  these  results  shows  that 
the  effect  of  including  normal  pressure  gradients  is  relatively  minor  and 
has  less  of  an  impact  on  the  results  than  the  change  in  the  turbulence 
model  discussed  previously.  A  similar  conclusion  is  reached  by  examining 
the  velocity  profiles.  Typical  results  are  shoim  in  Fig.  30  for  the  flow 
just  upstream  and  downstream  of  the  shock  wave  in  which  only  minor  changes 
occur  in  the  velocity  profile  when  the  normal  pressure  gradient  is  deleted 
from  the  analysis.  Although  the  present  treatment  of  normal  pressure 
gradients  is  an  approximate  one,  it  appears  from  these  calculations  that 
the  dominant  viscous  effect  in  the  prediction  of  separated  transonic  shock- 
wave,  boundary- layer  interaction  is  the  inclusion  of  displacement  thickness 
interaction.  The  solutions  shown  in  Figs.  29  and  30  demonstrate  that 
meaningful  results  can  be  obtained  with  IBLT  without  the  inclusion  of 
normal  pressure  gradients.  In  the  future  work  should  be  directed  toward 
additional  data  comparisons  over  a  range  of  transonic  shock  strengths  to 
further  substantiate  this  conclusion. 


CONCLUDING  REMARKS 


In  this  report  a  technique  based  on  a  two- layer  model  has  been 
presented  for  the  analysis  of  transonic  shock-wave,  boundary-layer  interac¬ 
tion  with  particular  emphasis  on  turbulent  separated  flow.  In  this  tech¬ 
nique,  which  is  referred  to  as  Interacting  Boundary  Layer  Theory  (IBLT),  a 
finite  difference  solution  of  the  viscous  layer  equations,  expressed  in  a 
detect  form,  is  iteratively  solved  with  the  finite  difference  form  of  the 
inviscid  stream  function-vorticity  equation.  Normal  pressure  gradients  and 
embedded  shock  effects  have  been  included  in  the  viscous  layer  analysis  by 
assuming  that  the  pressure  in  this  region  is  the  same  as  that  of  the 
inviscid  flow  above  the  displacement  thickness.  A  number  of  IBLT 
calculations  for  both  attached  and  separated  flow  along  with  detailed 
comparisons  with  experimental  data  lead  to  the  conclusions  presented  below. 
In  addition  the  rotational  inviscid  flow  model  of  shock,  boundary- layer 
interaction  has  been  implemented  through  the  use  of  the  present  inviscid 
stream  function-vorticity  method.  Conclusions  are  also  drawn  from 
calculations  made  with  this  inviscid  interaction  approach. 


The  major  conclusion  which  is  dra%ra  from  the  present  work  is  that 
displacement  thickness  interaction  is  the  dominant  effect  which  must  be 
included  in  a  viscous-inviscid  interaction  scheme  for  the  prediction  of 
transonic  shock-wave,  boundary- layer  interaction  where  the  shock  is  of 
sufficient  strength  to  separate  the  turbulent  boundary  layer.  Comparisons 
with  the  e  1 .4  case  measured  by  Kooi  indicate  that  for  this  separated 
flow,  the  computed  results  shows  greater  sensitivity  to  the  turbulence 
model  than  to  whether  or  not  normal  pressure  gradients  are  included  in  the 
analysis.  In  this  case  the  pressure  was  found  to  vary  only  by  about  10 
percent  across  the  viscous  layer;  hence  it  is  not  surprising  that  the  IBLT 
calculations  proved  to  be  relatively  insensitive  to  whether  or  not  this 
effect  was  included.  In  the  future  additional  data  comparisons  should  be 
made,  particularly  as  the  shock  strength  is  increased,  to  further  substan¬ 
tiate  this  conclusion.  Unfortunately,  most  of  the  data  for  transonic 
shock-wave,  boundary- layer  interaction  has  been  taken  in  small  wind  tunnels 
where  it  is  necessary,  as  was  found  in  the  present  analysis  of  Kooi's  data, 
to  include  the  effect  of  boundary  layer  blockage. 

Favorable  comparisons  have  been  obtained  with  the  present  IBLT  for 
Kooi's  data  for  the  pressure,  displacement  thickness,  and  velocity  profiles 
in  which  a  modified  Cebeci-Smith  turbulence  model  has  been  used.  This 
modification  which  is  a  straightforward  change  in  the  length  scale  in  the 
outer  region,  provides  encouragement  that  a  simple  eddy  viscosity  model  may 
be  adequate  for  the  prediction  of  many  of  the  flow  proo'^-ties  in  a 
turbulent  separated  flow.  However  the  skin  friction  ano  momentum  thickness 
comparisons  indicate  that  further  work  is  needed  to  improve  this  model  in 
the  inner  or  wall  region. 


The  present  calculations  demonstrate  that  IBLT  can  produce  meaningful 
overall  results  for  this  complex  flow  field.  However,  the  details  of  the 
lambda  shock  pattern,  generally  observed  when  the  flow  separates,  are  not 
present  in  the  IBLT  calculation.  Resolution  of  this  flow  structure  would 


39 


require  inclusion  of  the  exact  shock  jump  conditions  and  the  use  of  an 
extremely  fine  computational  mesh  in  order  to  delineate  the  details  of  the 
triple  shock  point  and  the  entropy  and  total  pressure  jump  which  occur 
across  the  shear  layer  emanating  from  this  point. 

Favorable  comparisons  which  have  been  made  with  the  asymptotic  theory 
of  Melnik  and  Grossman  indicate  that  the  present  transonic  stream  function- 
vorticity  approach  can  be  used  as  an  inviscid  rotational  flow  model  of 
shock  interaction  if  the  flow  remains  attached.  Calculations  are  also 
included  which  show  that  this  model,  as  expected,  is  completely  inadequate 
if  the  flow  separates.  In  this  case,  viscous  effects,  as  demonstrated  by 
the  sensitivity  of  the  present  results  to  the  displacement  thickness  and 
the  turbulence  model,  must  be  included  in  the  interaction  analysis. 

A  favorable  comparison  was  also  obtained  with  the  attached  flow 
results  of  Melnik  and  Grossman  with  the  present  IBLT,  without  the  inclusion 
of  normal  pressure  gradients.  The  advantage  of  this  approach  over  the  ro¬ 
tational  inviscid  flow  theory  is  that  the  entire  boundary  layer  is  solved 
at  once  without  breaking  it  into  numerous  layers  which  require  several 
solutions  to  obtain  the  entire  solution. 
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Fig.  1  Normal  shock  wave  impinging  on  a  turbulent  boundary  layer 
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Fig.  2  Viscous  —  inviscid  interaction  fiow  modeis 
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Fig.  6  Viscous-inviscid  matching  of  veiocity  profiie 
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Fig.  13  Comparison  of  wall  prassure  from  present  rotational  inviscid  and  interacting 
boundary  layer  theory  (IBLT)  with  asymptotic  theory  of  Melnik  and  Grossman 


Fig.  17  Ffm  air  iBLT  caicuiation  of  normai  shock*wava,  boundary-iayer  interaction 
a)  Wali  preasure 
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Fig.  18  UpstrMfn  profiles  for  inviscid  rotational  flow  analysis  of  normal  shock*wave, 
boundary-layer  interaction  at  Mgo  «1.4 
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Fig.  21  Comparison  of  results  from  IBLT  (py^tO)  with  experimental  data  for  transonic 
normal  shock'Wave,  boundary*layer  interaction 
a)  Wall  pressure 
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Fig.  21  Comparison  of  rosulta  from  IBLT  (Py^O)  with  axparimantal  data  for  transonic 
normal  shock-wava,  boundary-layar  intaraction 
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Fig.  21  Comparison  of  rasults  from  IBLT  (Py^O)  with  axparimantal  data  for  transonic 

normal  shock-wava,  boundary<layar  Intoraction 
a)  Shapa  factor 


Fig.  22  Comparison  of  veiocity  profiies  from  iBLT  (py^^O)  with  experimentai  data 
for  Mgp  =1.4,  Re 05,  =2.97  x  10® 
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Variation  of  pressure  across  boundary  iayer  from  equivaient  inviscid  flow  (EIF)  for 
transonic  normal  shock-wave,  boundary-layer  interaction 


PROFILE  LOCATIONS 


Velocity  from  real  viscous  flow  (RVF)  and  equivalent  inviscid  flow 
transonic  normal  shock-wave,  boundary-layer  interaction 
a)x=  -2.5,  2.0,  2.9 


PROFILE  LOCATIONS 


Velocity  from  real  viscous  flow  (RVF)  and  equivalent  inviscid  flow  (EIF) 
for  transonic  normal  shock-wave  boundary  layer  interaction 
b)x  =  3.5,  3.8,  and  4.4 


Fig.  26  Velocity  from  real  viscous  flow  (RVF)  and  equivalent  inviscid  flow  (EIF) 
for  transonic  normal  shock-wave,  boundary-layer  interaction 
c)  x  =  5.6,  8.0,  and  11.0 
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Stream  function  profiles  from  real  viscous  flow  (RVF)  and  equivalent  inviscid 
flow  (EIF)  for  transonic  normal  shock-w&ve,  boundary-layer  interaction 
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Fig.  29  Comparison  of  computed  resuits  from  iBLT  with  and  without  normal  pressure 
gradients  across  viscous  layer  for  transonic  normal  shock>wave,  boundary-layer 
interaction 
a)  Wall  pressure 


Fig.  29  Comparison  of  computed  resuits  from  IBLT  with  and  without  normai  pressure 
gradients  across  viscous  layer  for  transonic  normal  shock-wave,  boundary-layer 
interaction 

b)  Displacement  thickness 
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Fig.  30  Comparison  of  veiocity  profiies  from  e 
without  normai  pressure  gradients  acr 


CALCULATIONS  OF  SEPARATION  BUBBLES  USING 
BOUNDARY-LAYER-TYPE  EQUATIONS 


A.  Halim*  and  M.  Hafez** 


1.  INTRODUCTION 

Incompressible  laminar  flows  with  small  separation  bubbles 
have  been  successfully  simulated  using  Navier  Stokes  equations 
for  some  years  (see  Reached))*  For  example,  Briley d)  and 
Leal(3)  used  stream  function-vorticity  formulation  while  Ghia 
et  al(^)  used  primitive  variables.  Simplified  equations  of  the 
first  formulation  (where  the  streamwise  viscous  terms  are  ne¬ 
glected)  were  used  by  ^ia  and  Davisd)^  Werle  and  Bemsteind)^ 
Ghia,  Ghia  and  Teschd),  and  Inoued).  Similarly,  partially 
parabolized,  or  semi-elliptic  equations  in  tenas  of  primitive 
variables  were  used  by  Mahgoub  and  Bradshawd),  and  Rubin dO), 

In  the  first  approach,  careful  treatment  of  the  vorticity 
boundary  condition  at  the  wall  is  required,  while  conservation 
of  mass  is  usually  difficult  to  achieve  in  the  second  approach 
(see  Orszag  and  Israeli dl)). 

Recently,  boundary  layer  equations  have  been  used  to  cal¬ 
culate  separated  flows.  There  are  three  main  difficulties  with 
classical  boundary  layer  methods.  First,  it  is  well  knoxro  that 
when  the  pressure  is  prescribed,  the  equations  admit  a  singular 
solution  (see  the  reviews  by  Brown  and  Stewartsond2)  ^nd 
Williams d3))^  By  differentiating  the  x-momentum  equation 
twice  with  respect  to  y  and  using  the  continuity  equation,  the 
shear  stress  at  the  wall  is  given  by 
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Etence,  unless,  Uyyyy  vanishes  at  the  separation  point,  the 
shear  stress  has  a  square  root  singularity.  More  precisely, 
Goldstein showed  that  the  function  Lim  (Re^/^x)  which  is 

Re-Mo 

generally  boxinded  upstream,  is  nonanalytlc  at  the  point  of  zero 
skin  friction.  There  is  no  contradiction  between  this  theory 
and  Dean's  results which  show  that  the  Mavler-Stokes  equa¬ 
tions  are  analytic  there.  The  limiting  processes  are  different. 

The  second  difficulty  of  direct  boundary  layer  calcula¬ 
tions  is  the  possibility  of  nonunique  solutions.  This  is  clear 
from  the  similar  solution.  It  has  also  been  identified  for 
non-similar  solutions  (see  Murphy  and  Kingd^)  and  the  referen¬ 
ces  therein). 

Finally,  the  marching  procedure  is  not  stable  for  reversed 
flows.  Flugge  Lotz  and  Reyhner^^^)  neglected  the  streamwise  con¬ 
vective  term  when  u  is  negative.  More  accurate  solutions  can 
be  obtained  by  either  using  the  appropriate  upwind  difference 
schemes  and/or  integrating  the  equations  in  the  appropriate 
direction  locally.  The  stability  of  such  calculations  can  be 
improved  by  solving  artificial  time-dependent  equations  where 
the  diagonal  dominance  is  enhanced  due  to  the  time-like  term. 

The  numerical  breakdown  of  classical  boundary  layer  methods 
near  a  separation  point  has  been  demonstrated  by  many  authors 
(e.g.  Wele  and  Davis (i®),  Klienberg  and  Steger^^^^,  and 
Fletcher  and  Dancey(20)).  xhe  solution  becomes  singular  and  it 
is  difficult  to  proceed  the  calculations  downstream  of  the 
separation  point.  Following  Goldstein's  analysis,  it  can  be 
shown  that  near  the  separation  point  x^, 

6*(x)  ■  '5*(Xg)  -  +  ...  (2) 

dT  d6* 

So  that  if  is  infinite  there*  so  is  Therefore  regular¬ 

ity  can  be  ensured  by  requiring  either  T(x)  or  5*(x)  to  be  re¬ 
gular.  Catherall  and  Mangler^^^^  were  the  first  to  obtain 
regiilar  solution  by  specifying  the  displacement  thickness. 

Carter demonstrated  that  results  of  Inverse  boundary  layer 
calculations  are  in  good  agreement  with  Navier  Stokes  solution. 
Similar  success  was  reported  by  Cebecl,  Keller,  and  Williams (23) 
using  a  nonlinear  eigenfunction  formulation.  On  the  other 
hand,  Klienberg  and  Steger(29)  and  Horton(2^)  specified  the 
shear  stress  in  their  calculations  and  obtained  regular  solu¬ 
tions.  Indeed,  Horton  concluded  that  the  condition  uyyyy*0,  at 
the  separation  point  was  satisfied  up  to  the  computational 
accuracy.  It  should  be  mentioned  that  specifying  the  displace¬ 
ment  thickness  leads  to  faster  convergence  and  it  seems  it 
avoids  the  nonuniqueness  problem. 


The  appropriate  dlsplaceoent  thickness  or  wall  shear  must 
be  obtained  as  part  of  the  overall  problem,  from  the  inter¬ 
action  between  the  boundary  layer  and  the  inviscid  flow.  Cal¬ 
culation  of  laminar  separation  bubbles  that  Includes  this  vis¬ 
cous  Invlscld  Interaction  has  been  given  by  Carter  and 
WomornC^^).  On  the  other  hand.  It  is  argued  that  due  to  the 
interaction,  the  mainstream  is  always  able  to  adjust  itself  to 
prevent  the  singularity  from  developing.  Briley  and 
McDonald (27)calculated  Che  leading  edge  separated  flows 
using  a  direct  procedure,  where  the  unsteady  boundary  layer 
equations  are  repeatedly  solved  until  a  steady  state  solution 
is  obtained;  after  each  time  step,  the  prescribed  pressure  is 
updated  from  thin  airfoil  theory,  thereby  accounting  for  the 
displacement  interaction.  The  same  problem  has  been  also  sol¬ 
ved  by  Kwon  and  Fletcher (28) ^  Cebeci  and  Schimke(29)^  Vatsa 

and  Carter(^(^'.  Davis  and  Werle^^^^ (^2)  viewed  the  interactive 
boundary  layer  as  an  alternative  to  the  asymptotic  (triple 
deck)  theory  of  Stewartson  and  Messiter.  Recently,  Carter 
LeBalleur(2^) ,  Wigton^^S)^  and  Veldman^^^)  introduced  differ¬ 
ent  coupling  procedures  to  accelerate  the  convergence  of  the 
interaction.  For  some  special  cases,  inviscid  rotational  flow 
(see  Taulbee  and  Robertson^^'' )  or  just  boundary  layer  (see 
Pletcher(^®))  provide  reasonable  approximate  solutions.  In 
general,  neither  is  uniformly  valid  and  the  interaction  is  re¬ 
quired.  The  interaction  reflects  the  boundary  value  nature  of 
the  problem  tdilch  is  important  even  for  supersonic  flows  where 
(space-)  marching  procedures  can  be  used  for  both  the  inviscid 
and  the  boundary  layer  parts  separately  (see  Garvine^^^^). 

In  this  paper,  numerical  calculations  of  separated  flows 
are  presented  where  boundary  layer  equations  in  terms  of  stream 
function,  with  different  bomdary  conditions  are  integrated. 

For  a  retarded  flow  on  a  flat  plate,  results  are  compared  with 
Briley's  and  Carter's  solutions.  More  accurate  models  account¬ 
ing  for  the  pressure  variation  across  the  viscous  layer  are 
discussed.  In  many  cases, the  effect  of  such  a  variation  is  neg¬ 
ligible. 

For  viscous  inviscid  interaction  problems,  a  patching  pro¬ 
cedure  is  used  where  the  unknowns  along  a  vertical  line  in  both 
domains  are  solved  simultaneously  leading  to  strong  coupling. 

In  this  method,  the  grid  should  represent  both  the  viscous  and 
the  inviscid  scales,  but  the  equations  are  simpler  than  the  full 
Navler  Stokes  (biharmonlc)  equation. 

2.  A  MODEL  PROBLEM 
(2) 

Briley'  '  suggested  to  calculate  retarded  flows  on  a  flat 
plate  using  Navier  Stokes  equations  and  compared  the  results 
with  Howarth's  boundary  layer  series  solution  where  it  breaks 
down  at  the  separation  point  as  predicted  by  Goldstein  analy¬ 
sis.  Later,  Carter(22)  prescribed  the  displacement  thickness 
from  Briley's  solution  and  demonstrated  that  the  inverse 


p>. 


boundary  layer  solution  Is  In  excellent  agreement  with  Briley's 
results.  Recently  many  authors  solved  the  same  problem  - 
Anarante  and  Cheng and  Cebecl  and  Stewartson^^^)  used 
Interactive  boundary  layer  methods  while  Inoue(8)  used  para¬ 
bolized  Navler  Stokes  equations.  We  chose  this  model  problem 
to  test  present  formulations. 


3.  BOUNDARY  LAYER  EQUATIONS 


In  terms  of  the  tangential  and  the  normal  velocity  compo¬ 
nents,  the  continuity  and  the  tangential  momentum  equations 
read: 


U  +  V  = 

0 

i 

X  y 

UU  +  vu 

X  V 

■  -  p  + 

*^x 

yy 


(3) 

(4) 


with  the  boundary  conditions  u*v=0  at  y=0,  and  u*Ug  at 
(edge  of  boundary  layer),  where  It  Is  assumed  that  px  does  not 
vary  with  y.  Equations  (3)  and  (4)  can  be  written  In  terms  of 
stream  function  and  vortlclty,  namely: 


Uu  +  V(^ 

X  y 


VO) 


yy 


-OJ 


yy 


(5) 

(6) 


with  the  boundary  conditions  ip^ipy^O  at  y“0  and  ij>y*Ug  and  a)=0  at 

y*ye- 


Another  form  of  the  above  equations  Is  simply: 


u  U/  +  V  \l/ 

yyx  yyy 


yyyy 

The  advantage  of  form  (7)  Is  discussed  next 


(7) 


First,  the  discretization  Is  straightforward;  (for  uni¬ 
form  mesh) : 


(8-1) 


(8-2) 


if  u>0  (8-3) 
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(8-4) 


if  u<0 

u  and  V  are  expressed  in  terms  of  tp: 


or 


(8-5) 

''c-D  “  ■  <*i+l,j  “  *1-1,3'^^^ 

(8-6) 

"b-d  “  -  '*1.3  -  *1-1.3>' 

(8-7) 

At  each  x“constant  line  the  unknowns  form  a  five  diagonal 
scalar  system  of  equations.  Both  wall  boundary  conditions  can 
be  easily  Incorporated  implicitly.  Similarly,  the  edge  bound¬ 
ary  conditions  can  be  Implemented  by  augmenting  the  system  with 
the  extra  relations , for  example  '{>y“Ug  and  li^yy^O.  In  this  case, 
the  non-zero  coefficients  of  the  matrix  are  shown  in  sketch  (1). 
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Sketch  (1)  -  Nonzero  Elements  of  the  System  Matrix 

The  resulting  equations  can  be  solved  efficiently  by  a  Gaussian 
elimination  procedure  especiallzed  to  five  diagonal  scalar  sys¬ 
tems  (an  extension  of  Thomas  algorithm).  For  attached  flows, 
assuming  an  initial  profile  is  given  and  with  a  backward  differ¬ 
ence  approximation  of  v,  equation  (7)  is  integrated  by  marching 
in  the  x-dlrectlon.  Few  iterations  are  required  at  each  sta¬ 
tion  to  account  for  the  nonlinearity.  For  separated  flows, 
centered  differences  are  used  for  v  and  a  downstream  boundary 
condition  seems  to  be  required.  Two  possibilities  are  tested, 
at  the  last  station,  a  backward  difference  approximation  is 
used  for  v  or  v  is  set  equal  to  zero.  Also,  global  (sweeps) 
rather  than  local  iterations  are  employed. 


A.  PARTIALLY  PARABOLIZED  NAVIER  STOKES  EQUATIONS 

In  this  model,  only  the  streamwise  viscous  terms  in  the 
full  equations  are  neglected,  namely. 


U  +  V  * 
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It  is  argued  that  the  term  Vyy  in  equation  (11)  is  also  of 
higher  order. 

Equations  (9-11)  can  be  written  in  terms  of  stream  func¬ 
tion  and  vorticity: 

tjoi  +  VO)  =  voj  (12) 

X  y  yy 

.  -  (jj  =4;  +4/  (13) 

XX  yy 

The  difference  between  equations  (12)  and  (13)  and  the  boundary 
layer  equations  (5)  and  (6)  is  the  term  li'xx  equation  (13). 
Inoue^®)  solved  equations  (12)- (13)  for  Briley's  problem  and 
obtained  close  results. 

Substituting  equation  (13)  into  equation  (12)  leads  to  a 
fourth  order  equation  in  ’4;; 

u(4'  +'+')+  v(4  +4^  )  =  v(i|'  +4'  )  (lA) 

'  XX  yy^X  '  XX  yy-'y  '■’^xx  TJ  Tf 

Centered  differences  are  always  used  for  'J'xx'*'’^yy »  differ¬ 

ences  are  used  for  the  streamwise  convection  operator  uSx-  The 
wall  boundary  conditions  are  4'*'i'y“0  and  at  the  outer  edge 
4^xx''^^yy*^  and  4^y*Ug. 

To  solve  equation  (14),  two  artificial  time  dependent  terms 
are  added,  04/^  and  B(4'xx'^'^^yy)t  enhance  diagonal  dominance. 
Again  a  five  diagonal  scalar  system  at  each  station  is  solved 
in  the  same  manner  the  boundary  layer  equation  (7) is  solved. 

At  the  downstream  boundary,  the  4'xx  equation  (14)  is 

neglected  and  the  equation  reduces  to  equation  (7).  Global 
iteration  is  always  needed. 

The  extra  terms  in  equation  (14)  and  not  in  equation  (7) 
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These  terns  are  equal  to  v^lch  has  been  neglected  in  boun¬ 
dary  layer  calculations.  For  many  cases,  a  boundary  layer  type 
equation  can  be  used,  where  the  term  p^^y  Is  lagged.  For  exam¬ 
ple,  the  pressure  can  be  obtained  from  a  triadiagonal  system 
resulting  from  the  discretization  of: 

'yy  ■  ■  V'y  (15) 

with  the  boundary  conditions: 

Py  ■  0  at  the  wall 

and  p  -  H  -  (u^  +  v^)  (16) 

at  the  outer  edge,  where  H  Is  a  constant.  In  equation  (13) , 
the  viscous  term  '^^yyy  has  been  neglected.  After  the  pressure 
is  calculated,  the  term  p^  Is  updated.  To  ensure  stability  of 
calculations,  the  pressure  Is  underrelaxed. 

5.  VISCOUS /INVISCID  INTERACTION 

The  viscous  terms  are  negligible  far  from  the  wall,  and 
the  flow  becomes  Invlscld.  There  are  several  methods  for  boun¬ 
dary  layer  and  Invlscld  calculations.  Four  possibilities  are 
considered : 

1.  using  Hilbert  integral  for  invlscld  flow  and  Integral 
boundary  layer  calculations, 

2.  using  Hilbert  Integral  and  finite  difference  boundary 
layer  calculations, 

3.  using  finite  differences  for  invlscld  flows  and  inte¬ 
gral  boundary  layer  calculations, 

4.  using  finite  differences  for  both  invlscld  and  boun¬ 
dary  layer  flows. 

In  our  opinion,  except  for  class  4,  the  Invlscld  and  viscous 
equations  should  be  solved  simultaneously  (see  for  example 
references  (42),  (36),  and  (43)).  A  coupling  procedure  Is 
needed  only  for  class  4,  where  at  each  step,  the  Invlscld  and 
the  boundary  layer  equations  are  solved  separately;  the  coupling 
procedure  update  the  boundary  conditions  for  both  calculations 
and  the  process  Is  repeated  until  convergence. 

An  alternative  formulation,  based  on  a  patching  (or  zonal) 
procedure  Is  tested  here.  One  requirement  for  the  present 
method  Is  the  existence  of  an  overlap  region,  where  the  equa¬ 
tions  on  both  sides  are  valid.  In  the  viscous  region,  Che  par¬ 
tially  parabolized  Navler  Stokes  equations  are  used  where  the 
vorticity  Is  convected  and  diffused.  On  the  other  hand,  the 


equation  for  the  Invlscld  region  is 


'I' _  +  It. 
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where 


(17) 


cj  •  (18) 

i. e. ,  the  vorticity  is  only  convected.  Hence,  the  patching 
line  should  be  in  a  region,  where  the  viscous  terns  are  negli¬ 
gible.  If  the  patching  line  is  chosen  far  from  the 
wall,  the  vorticity  there  vanishes  and  the  flow  is  Irrotatlonal. 
It  is  clear  that  this  is  a  special  case  and  the  Irrotationallty 
is  not  necessarily  required  in  the  present  method.  Numerical 
Impl&Dentation  is  simple.  At  each  vertical  line,  the  viscous 
terms  are  tested  (for  all  the  points),  when  they  are  relative¬ 
ly  small,  the  five  diagonal  system  is  replaced  by  a  tridiagonal 
system  resulting  from  equation  (17).  It  would  be  more  effici¬ 
ent,  to  have  a  dynamic  grid  generation  method,  where  the  grid 
points  are  distributed  according  to  the  proper  scales  based  on 
order  of  magnitude  analysis. 

If  one  wishes,  for  some  reason  (say  to  take  full  advantage 
of  the  different  scales  Involved),  to  solve  the  viscous  and  the 
Invlscld  problem  separately,  a  coupling  procedure  is  needed. 

Two  methods  have  been  tested,  the  first  is  a  direct  local 
coupling  and  the  second  is  a  semi-inverse  implicit  procedure  as 
shown  in  Sketches  (2)  and  (3)  respectively. 


Invlscld 


Viscous 


Sketch  (2)  -  A  direct  local  coupling  procedure 
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Viscous 
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Formula  ] 

Sketch  (3-a)  -  A  semi-inverse  implicit 
coupling  procedure. 
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unknowns  in  Inviscid  Region,  »  matrix  of  inviscid 
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unknowns  in  Viscous  Region,  ■  matrix  of  viscous  equa¬ 
tions 

unknowns  along  the  patching  line,  ■  matrix  of  coupling 
equations 


Sketch  (3-b)  -  Partitioning  of  the  total 
matrix  for  the  Implicit  coupling  pro¬ 
cedure 
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In  the  first  method,  the  patching  line  is  chosen  between  grid 
points  say  P  and  P+1.  from  the  viscous  region  is  used  as 

boundary  condition  for  the  Invlscid  calculations  which  provide 
and  4'p.f2  ^  boundary  conditions  for  the  viscous  calcula¬ 
tions.  Convergence  is  achieved  when  ipp,  and  V'p^.o 

same  (within  some  tolerance)  in  both  calculations.  The  second 
method  starts  with  an  initial  guess  for  p.  The  invlscid  prob¬ 
lem  is  solved  first,  followed  by  the  viscous  problem  using  the 
same  values  for  4^  as  well  as  'I'p^.^  from  the  invlscid  output. 

Next  4'p  is  updated  through  a  coupling  formula,  equation  (19), 
and  the  process  is  repeated.  The  coupling  formula  is  consistent 
with  the  invlscid  equation,  namely. 


,,n+l  o,n+l  ,  ,n+l  x#,  2 
^'^i+l,P  "  ^'*'i,P  ’'^'i-l,P^^^* 


+  “  ■  w(i,P) 


where  a  is  a  relaxation  parameter,  p^2.  is  obtained  from  the 
invlscid  calculations,  while  the  viscous  calculations  provide 
4(i^p_l.  A  trldlagonal  solver  is  used  to  obtain  If  the 

x-term  and  ui  are  neglected,  equation  (19)  reduces  to: 


+ 


(u,  -  u  ) 
i  V 


Equation  (20)  is  similar  to  Carter's  formula It  is  found, 
however,  that  equation  (19)  leads  to  faster  convergence. 


(As  a  side  remark,  we  notice  that,  for  transonic  calculations, 
equation  (19)  should  reflect,  the  local  character  of  the  flow, 
whether  subsonic  or  supersonic  and  in  this  case,  equation  (19) 
is  replaced  by: 
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where  6il>  is  4)  -  4^  and  6  and  a  are  relaxation  parameters  and 


d  is  an  artificial  density  to  produce  dissipation  in  the  super¬ 
sonic  region.  Transonic  calculations  will  be  reported  else- 
%fhere. ) 


6.  NUMERICAL  RESULTS 


Three  uniform  meshes  (in  x  and  y)  are  used:  Coarse  (35, 
77),  intermediate  (76,77),  and  fine  (151,77).  Howarth  profile 
is  specified  at  the  upstream  boundary.  Different  edge  boundary 
conditions  are  tested.  Five  cases  are  discussed  here. 


Case  1:  4'g*Ug(yg-5*)  t  and  4)^*0,  where  6*  is  prescribed  from 


Brilev's  solution. 


• 


Case  2:  -u  (y  -6*),  and  ip  “u  . 

e  e'-^e  ’  ye 


Case  3:  "u  ,  and  \1)  “0,  where  u  is  the  output  of  Case  1. 

y  *  yy  ^ 

Case  4:  i|f  "u  ,  and  \p  >0. 

ye’  ^yy 

The  results  shown  for  Cases  1  to  4  are  for  intermediate  mesh. 

Case  5:  Same  as  Case  4  but  using  fine  mesh. 

Case  1  is  the  same  as  Carter’ s  inverse  boundary  layer  calcula¬ 
tions.  Case  2  has  a  hybrid  boundary  condition  where  both  ij/  and 
are  prescribed.  As  expected,  the  solution  is  regular  since 
by  specifying  il>,  v  is  guaranteed  to  be  bounded.  We  notice  that 
if  and  'py  are  prescribed  at  the  wall  and  at  y=yet  an  approxi¬ 
mate  solution,  with  closed  streamlines  is  obtained  using  a  cubic 
polynomial  ('i'yyyy“0).  Case  3  is  interesting  since  it  demon¬ 
strates  that  the  calculations  are  stable  and  the  solution  is 
unique  and  regular.  The  results  of  Cases  4  and  5  are  somewhat 
surprising  since  this  is  similar  to  the  classical  direct  boun¬ 
dary  layer  calculations.  There  are,  however,  some  differences, 
since  centered  differences  are  used  for  v  and  a  downstream 
boundary  condition  is  required.  These  results  are,  of  course, 
for  finite  Reynolds  number  and  finite  grid  size  (the  scheme  is 
first  order  accurate).  Indeed,  if  backward  differences  are 
used  for  V,  the  solution  becomes  singular  as  approaching  the 
separation  point  and  the  calculation  is  not  stable  downstream 
of  the  zero  skin  friction  point. 

Figure  (1)  shows  the  streamlines  calculated  by  Briley  and 
by  Carter.  In  figure  (2),  the  corresponding  streamlines  of  the 
five  cases  discussed  above  are  presented.  The  rates  of  conver¬ 
gence  are  plotted  in  figure  (3).  We  notice  that  the  convergence 
is  faster  when  ip  is  prescribed.  Comparing  Case  1  to  Carter's 
calculation,  it  is  at  least  twice  as  fast.  Figures  (4)  and  (5) 
show  the  displacement  thickness  and  the  wall  shear  distributions. 

The  discrepancy  between  the  boundary  layer  solutions  of 
Cases  (4)  and  (5)  and  the  Navler  Stokes  solution  may  be  due  to 
the  different  boundary  conditions.  In  Briley’s  calculation, 
<(<xx+’('yy*“‘*’"0  is  used,  while  the  boundary  layer  approximation  of 
the  vortlclty,  \pyy  is  set  equal  to  zero  in  Cases  4  and  5.  From 
the  streamlines  shown  in  figure  (1),  it  is  clear  that  at  y*ye> 
ipn  does  not  vanish. 

If  viscous  inviscld  Interaction  is  allowed,  the  results  do 
not  agree  with  Briley’s,  since  in  his  study,  it  was  assumed  chat 
"along  the  outer  flow  boundary  u  was  available  from  an  inviscid 
solution.  In  the  resulting  Navier  Stokes  solution,  the  v- 
component  of  velocity  along  the  outer-flow  boundary  generally 
will  not  agree  with  that  from  the  inviscld  solution  used  to 
prescribe  the  distribution  of  u  along  that  boundary,  and  the 
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mismatch  in  v  Is  an  indication  of  interaction  between  the  In- 
viscid  and  viscous  solutions  which  has  been  neglected.  To 
account  for  the  interaction,  it  would  be  necessary  to  success¬ 
ively  recompute  the  Invlscld  and  Navier-Stokes  solutions, 
siloing  in  some  manner  for  the  Influence  of  one  upon  the 
other,  until  the  two  solutions  no  longer  change  significantly." 

The  present  viscous  Invlscld  Interaction  method  described 
above,  does  not  produce  any  separation  for  the  cases  Briley 
solved.  In  a  recent  paper,  by  Cebeci  and  Stewartson^^^) ,  it 
was  shown  that,  based  on  interactive  boundary  layer  calcula¬ 
tions,  the  critical  value  of  Xg  (a  parameter  defining  the 
comer  point  of  the  external  velocity)  to  induce  separation  is 
0.215  for  Re>10^/A8.  Comparable  results,  for  the  displacement 
thickness  are  shown  for  XqO.IS  (corresponding  to  Briley's 
solution  for  Xo*0.202).  They  found  that  there  is  a  maximum 
value  of  Xo,  after  which  the  numerical  procedure  breaks  down. 
Furthermore,  this  maximum  value  Is  a  decreasing  function  of 
Reynolds  number  and  seems  to  approach  the  value  0.12  predicted 
by  Howarth.  It  was  concluded  that  there  Is  a  bound  to  the  use¬ 
fulness  of  interactive  boundary  layer  theory,  and  once  It  Is 
exceeded  the  theory  In  some  sense  goes  sour;  possibly  the  glo¬ 
bal  flow  properties  then  rapidly  change  over  to  those  corres¬ 
ponding  to  Klrchoff  free-streamline  flow.* 

Cebeci  and  Stewartson  used  Hilbert  Integral  for  Invlscld 
flow  and  Keller's  box  scheme  for  the  boundary  layer  calcula¬ 
tion.  Their  results  are  compared  to  the  present  viscous 
Invlscld  calculations  where  finite  differences  are  used  for  both 
regions,  as  shown  in  figure  (6).  The  agreement  is  reasonable, 
taking  Into  consideration  the  finite  domain  of  integration  of 
the  finite  difference  calculations  (130x115  points  are  used 
with  Xg-1. 2,  and  ye*^  ■  8.11).  The  performance  of  the  pre¬ 
sent  method  for  the  cases  where  the  interactive  boundary  layer 
calculations  diverge  is  the  subject  of  further  investigation. 

7.  CONCLUDING  REMARKS 

The  boundary  layer  equations.  In  terms  of  a  stream  func¬ 
tion  (differentiated  to  yield  a  fourth  order  equation)  is 
Integrated,  marching  in  the  main  stream  direction.  A  special 
Gaussian  elimination  procedure  is  applied  to  a  scalar  five  dla:- 
gonal  system  of  equations  at  each  step,  which,  in  a  sense.  Is 
similar  to  Davis  Coupled  Method^^^)  (where  continuity  and  mom¬ 
entum  equations  are  solved  together).  The  present  scheme  is 
not,  however,  as  contpact  as  Cranck-Nicholson' s  or  Keller's  box 
schemes (^5).  Nevertheless,  for  boundary  layer  type  calcula¬ 
tions,  a  smooth  stretching  transformation  is  usually  employed 
and  accuracy  Is  not  sacrificed  for  the  resulting  nonuniform 
mesh. 


Solutions  of  retarded  flows  are  presented  for  different 
outer  boundary  conditions,  namely  direct.  Inverse  and  hybrid. 

*Navler  Stokes  solutions  similar  to  those  of  Briley  (but  with 
higher  Re.  x^,  and  y  )  were  successfully  calculated  by 
Murphy * 
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Regular  solutions  at  the  separation  and  reattachment  points  are 
obtained.  Since  centered  differences  are  used  for  ^ 

downstream  boundary  condition  is  required  and  thus  nonunique¬ 
ness  as  well  as  singularities  are  avoided,  at  least  for  this 
formulation  with  the  grid  size  and  the  Reynolds  numbers  con¬ 
sidered. 

To  allow  for  viscous  Invlscld  interaction,  a  patching  pro¬ 
cedure  is  tested  inhere  a  partially  parabolized  Navler-Stokes 
equation  for  is  solved  simultaneously  with  the  invlscld 
Poisson  equation  via  a  vertical  line  relaxation.  For  each 
line,  a  coupled  system  (trldlagonal  in  the  invlscld  region  and 
five  diagonal  in  the  viscous  region)  is  solved  by  efficient 
Gaussian  elimination.  Artificial  time-dependent  terms  are 
added  to  enhance  the  stability  of  calculations. 

It  seems  that  the  present  stream  function  formulation 
leading  to  a  single  equation  Is  advantageous  for  discretization 
and  convergence.  In  ^is  work,  only  small  separation  bubbles 
are  calculated.  Simulation  of  large  regions  of  massive  separa¬ 
tion  would  be  of  great  Interest;  for  such  cases,  however,  the 
flow  In  reality  may  not  remain  laminar. 
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Figure  (5)  Vail  Shear  distribution  of  present  boundary 
layer  calculations 


2  .2  .6  .6  l.C  l.Z 


x/L 

Figure  (6)  Displacement  thickness  distribution  of  present 
viscous  inviscid  interactive  calculations 


